Fitting a Spatial Factor Multi-Species Occupancy Model

Data from: Plantas-Putumayo, Tamshiyacu-Tahuayo, Yasuni, Madidi

model
code
analysis
453 sites, 24 mammal species
Authors

German Forero

Robert Wallace

Galo Zapara-Rios

Emiliana Isasi-Catalá

Diego J. Lizcano

Published

August 16, 2025

Single-species occupancy models

Single-species occupancy models (MacKenzie et al., 2002) are widely used in ecology to estimate and predict the spatial distribution of a species from data collected in presence–absence repeated surveys. The popularity of these models stems from their ability to estimate the effects of environmental or management covariates on species occurrences, while accounting for false-negative errors in detection, which are common in surveys of natural populations.

Multi-species occupancy models

Multi-species occupancy models for estimating the effects of environmental covariates on several species occurrences while accounting for imperfect detection were build on the principle of the Single-species occupancy model, and developed more than 20 years ago.

Multi-species occupancy models are a powerful tool for combining information from multiple species to estimate both individual and community-level responses to management and environmental variables.

Objective

We want to asses the effect of protected areas on the occupancy of the species. We hypothesize that several mammal species have benefited of the conservation actions provided by protected areas, so we expect a decreases in occupancy from the interior of the protected area to the exterior as a function of the distance to the protected area border.

To account for several ecological factors we also tested and included elevation and forest integrity index as occupancy covariates.

We used the newly developed R package spOccupancy as approach to modelling the probability of occurrence of each species as a function of fixed effects of measured environmental covariates and random effects designed to account for unobserved sources of spatial dependence (spatial autocorrelation).

Code
library(grateful) # Facilitate Citation of R Packages
library(readxl) # Read Excel Files
library(DT) # A Wrapper of the JavaScript Library 'DataTables'
library(sf) # Simple Features for R
library(mapview) # Interactive Viewing of Spatial Data in R
library(maps) # Draw Geographical Maps
library(tmap) # Thematic Maps
library(terra) # Spatial Data Analysis
library(elevatr) # Access Elevation Data from Various APIs

# library(rjags) # Bayesian Graphical Models using MCMC 
library(bayesplot) # Plotting for Bayesian Models # Plotting for Bayesian Models
library(tictoc) # Functions for Timing R Scripts, as Well as Implementations of "Stack" and "StackList" Structures 
library(MCMCvis) # Tools to Visualize, Manipulate, and Summarize MCMC Output
library(coda) # Output Analysis and Diagnostics for MCMC
library(beepr) # Easily Play Notification Sounds on any Platform 
library(snowfall) # Easier Cluster Computing (Based on 'snow')

#library(ggmcmc)
library(camtrapR) # Camera Trap Data Management and Preparation 
library(spOccupancy) # Single-Species, Multi-Species, and Integrated Spatial Occupancy
library(tidyverse) # Easily Install and Load the 'Tidyverse'

Fitting a Multi-Species Spatial Occupancy Model

We use a more computationally efficient approach for fitting spatial multi-species occupancy models. This alternative approach is called a “spatial factor multi-species occupancy model”, and is described in depth in Doser, Finley, and Banerjee (2023). This newer approach also accounts for residual species correlations (i.e., it is a joint species distribution model with imperfect detection). The simulation results from Doser, Finley, and Banerjee (2023) show that this new alternative approach outperforms, or performs equally to spMsPGOcc(), while being substantially faster.

The Latent factor multi-species occupancy model is described in detail here

Data from Plantas Medicinales Putumayo

Here we use the table: - COL-18-Putumayo2023

Code
source("C:/CodigoR/WCS-CameraTrap/R/organiza_datos_v3.R")

AP_PlantasMed <- read_sf("C:/CodigoR/Occu_APs/shp/PlantasMedicinales/WDPA_WDOECM_May2025_Public_555511938_shp-polygons.shp")

### Ecu 17, Ecu 18, ECU 20, Ecu 22  

# load data and make array_locID column
Col_18 <- loadproject("F:/WCS-CameraTrap/data/BDcorregidas/Colombia/COL-018-Putumayo2023_WCS_WI.xlsx") |> mutate(array_locID=paste("Col_18", locationID, sep="_"))
60 cameras in Cameras. 
 60 cameras in Deployment. 
 60 deployments in Deployment. 
 60 points in Deployment. 
 51 cameras in Images. 
 51 points in Images. 
[1] "dates ok"
year: 2024 
 Jaguar_Design: no year: 2023 
 Jaguar_Design: no 
Code
# get sites
Col_18_sites <- get.sites("F:/WCS-CameraTrap/data/BDcorregidas/Colombia/COL-018-Putumayo2023_WCS_WI.xlsx")



# get elevation map
elevation_18 <- rast(get_elev_raster(Col_18_sites, z = 10)) #z =1-14
# bb <-  st_as_sfc(st_bbox(elevation_17)) # make bounding box 




# extract covs using points and add to _sites
covs_Col_18_sites <- cbind(Col_18_sites, terra::extract(elevation_18, Col_18_sites))
# covs_Col_17_sites <- cbind(Col_17_sites, terra::extract(elevation_17, Col_17_sites))


# get which are in and out
covs_Col_18_sites$in_AP = st_intersects(covs_Col_18_sites, AP_PlantasMed, sparse = FALSE)
# covs_Col_17_sites$in_AP = st_intersects(covs_Col_17_sites, AP_Yasuni, sparse = FALSE)




# make a map
mapview (elevation_18, alpha=0.5) + 
  mapview (AP_PlantasMed, color = "green", col.regions = "green", alpha = 0.5) +
  mapview (covs_Col_18_sites, zcol = "in_AP", col.regions =c("red","blue"), burst = TRUE) 
Code
  # mapview (covs_Col_17_sites, zcol = "in_AP", col.regions =c("red","blue"), burst = TRUE) +

Data fom Bajo Madidi and Heath

Here we use the tables: - BOL-008a - BOL-008b - BOL14a - BOL14b

Code
source("C:/CodigoR/WCS-CameraTrap/R/organiza_datos_v3.R")

Area_Madidi <- read_sf("C:/CodigoR/Occu_APs/shp/Area_Madidi/WDPA_WDOECM_Jul2025_Public_303894_shp-polygons.shp")

Madidi_NP <- read_sf("C:/CodigoR/Occu_APs/shp/Madidi_NP/WDPA_WDOECM_Jul2025_Public_98183_shp-polygons.shp")
#AP_Tahuayo <- read_sf("C:/CodigoR/Occu_APs/shp/Tahuayo/WDPA_WDOECM_Jun2025_Public_555555621_shp-polygons.shp")

# load data and make array_locID column
#Bol_Pacaya <- read_excel("F:/WCS-CameraTrap/data/BDcorregidas/Peru/PER-003_BD_ACRCTT_T0.xlsx", sheet = "Image")


Bol_Madidi_1 <- loadproject("F:/WCS-CameraTrap/data/BDcorregidas/Bolivia/Guido/BOL-008a.xlsx") |> mutate(Point.x=as.character(paste("M1",Point.x, sep = "-"))) |> mutate(Point.y=as.character(paste("M1",Point.y, sep = "-")))
65 cameras in Cameras. 
 61 cameras in Deployment. 
 61 deployments in Deployment. 
 31 points in Deployment. 
 61 cameras in Images. 
 31 points in Images. 
[1] "dates ok"
year: 2005 
 Jaguar_Design: yes 
Code
Bol_Madidi_2 <- loadproject("F:/WCS-CameraTrap/data/BDcorregidas/Bolivia/Guido/BOL-008b.xlsx") |> mutate(Point.x=as.character(paste("M2",Point.x, sep = "-"))) |> mutate(Point.y=as.character(paste("M2",Point.y, sep = "-")))
60 cameras in Cameras. 
 52 cameras in Deployment. 
 52 deployments in Deployment. 
 30 points in Deployment. 
 51 cameras in Images. 
 30 points in Images. 
[1] "dates ok"
year: 2005 
 Jaguar_Design: yes 
Code
Bol_Madidi_3 <- loadproject("F:/WCS-CameraTrap/data/BDcorregidas/Bolivia/Guido/BOL-014a.xlsx")|> mutate(Point.x=as.character(paste("M3",Point.x, sep = "-"))) |> mutate(Point.y=as.character(paste("M3",Point.y, sep = "-")))
47 cameras in Cameras. 
 44 cameras in Deployment. 
 44 deployments in Deployment. 
 26 points in Deployment. 
 44 cameras in Images. 
 26 points in Images. 
[1] "dates ok"
year: 2009 
 Jaguar_Design: yes 
Code
Bol_Madidi_4 <- loadproject("F:/WCS-CameraTrap/data/BDcorregidas/Bolivia/Guido/BOL-014b.xlsx")|> mutate(Point.x=as.character(paste("M4",Point.x, sep = "-"))) |> mutate(Point.y=as.character(paste("M1",Point.y, sep = "-")))
29 cameras in Cameras. 
 22 cameras in Deployment. 
 22 deployments in Deployment. 
 15 points in Deployment. 
 22 cameras in Images. 
 15 points in Images. 
[1] "dates ok"
year: 2009 
 Jaguar_Design: yes 
Code
# get sites
Bol_Madidi_sites1 <-  Bol_Madidi_1 |> 
   dplyr::bind_rows(Bol_Madidi_2) |> 
   dplyr::bind_rows(Bol_Madidi_3) |> 
   dplyr::bind_rows(Bol_Madidi_4) |> 
  select("Latitude", "Longitude", "Point.x" 
) |> dplyr::distinct( )  

Bol_Madidi_sites <- sf::st_as_sf(Bol_Madidi_sites1, coords = c("Longitude","Latitude"))
st_crs(Bol_Madidi_sites) <- 4326


# get elevation map
elevation_17 <- rast(get_elev_raster(Bol_Madidi_sites, z = 9)) #z =1-14
# bb <-  st_as_sfc(st_bbox(elevation_17)) # make bounding box 



# extract covs using points and add to _sites
covs_Bol_Madidi_sites <- cbind(Bol_Madidi_sites, terra::extract(elevation_17, Bol_Madidi_sites))
# covs_Ecu_17_sites <- cbind(Ecu_17_sites, terra::extract(elevation_17, Ecu_17_sites))


# get which are in and out
covs_Bol_Madidi_sites$in_AP = st_intersects(covs_Bol_Madidi_sites, Madidi_NP, sparse = FALSE)
# covs_Ecu_17_sites$in_AP = st_intersects(covs_Ecu_17_sites, AP_Machalilla, sparse = FALSE)


# make a map
mapview (elevation_17, alpha=0.7) + 
  mapview (Madidi_NP, color = "green", col.regions = "green", alpha = 0.5) +
  mapview (Area_Madidi, color = "yellow", col.regions = "yellow", alpha = 0.5) +
  #mapview (AP_Tahuayo, color = "green", col.regions = "green", alpha = 0.5) +
  mapview (covs_Bol_Madidi_sites, zcol = "in_AP", col.regions =c("red","blue"), burst = TRUE) 

Tamshiyacu Tahuayo data

Here we use the tables: - PER-003_BD_ACRCTT_T0.xlsx - PER-002_BD_ACRTT-PILOTO.xlsx

Code
source("C:/CodigoR/WCS-CameraTrap/R/organiza_datos_v3.R")

AP_Pacaya <- read_sf("C:/CodigoR/Occu_APs/shp/PacayaSamiria/WDPA_WDOECM_Jun2025_Public_249_shp-polygons.shp")

AP_Tahuayo <- read_sf("C:/CodigoR/Occu_APs/shp/Tahuayo/WDPA_WDOECM_Jun2025_Public_555555621_shp-polygons.shp")

# load data and make array_locID column
Per_Tahuayo_piloto <- loadproject("F:/WCS-CameraTrap/data/BDcorregidas/Peru/PER-002_BD_ACRTT-PILOTO.xlsx")
50 cameras in Cameras. 
 50 cameras in Deployment. 
 50 deployments in Deployment. 
 50 points in Deployment. 
 50 cameras in Images. 
 50 points in Images. 
[1] "dates ok"
year: 2015 
 Jaguar_Design: no 
Code
Per_Tahuayo <- loadproject("F:/WCS-CameraTrap/data/BDcorregidas/Peru/PER-003_BD_ACRCTT_T0.xlsx")
85 cameras in Cameras. 
 84 cameras in Deployment. 
 84 deployments in Deployment. 
 84 points in Deployment. 
 84 cameras in Images. 
 84 points in Images. 
[1] "dates ok"
year: 2016 
 Jaguar_Design: no 
Code
FLII2016 <- rast("C:/CodigoR/WCS_2024/FLI/raster/FLII_final/FLII_2016.tif")

# get sites
# Per_Pacaya_sites <- get.sites("F:/WCS-CameraTrap/data/BDcorregidas/Peru/PER-003_BD_ACRCTT_T0.xlsx")

# get sites
Per_Tahuayo_sites1 <- Per_Tahuayo |> select("Latitude",
                                    "Longitude",
                                    "Camera_Id" 
                                    ) |> dplyr::distinct( )  

Per_Tahuayo_sites <- sf::st_as_sf(Per_Tahuayo_sites1, coords = c("Longitude","Latitude"))
st_crs(Per_Tahuayo_sites) <- 4326



# get Pacaya sites
Per_Tahuayo_piloto_sites1 <- Per_Tahuayo_piloto |> select("Latitude",
                                    "Longitude",
                                    "Camera_Id" 
                                    ) |> dplyr::distinct( )  

Per_Tahuayo_piloto_sites <- sf::st_as_sf(Per_Tahuayo_piloto_sites1, coords = c("Longitude","Latitude"))
st_crs(Per_Tahuayo_piloto_sites) <- 4326




# get elevation map
elevation_PE <- rast(get_elev_raster(Per_Tahuayo_sites, z = 10)) #z =1-14
# bb <-  st_as_sfc(st_bbox(elevation_17)) # make bounding box 





# extract covs using points and add to _sites
covs_Per_Tahuayo_sites <- cbind(Per_Tahuayo_sites,
                                terra::extract(elevation_PE,
                                               Per_Tahuayo_sites)
                                )

# extract covs using points and add to _sites
covs_Per_Tahuayo_piloto_sites <- cbind(Per_Tahuayo_piloto_sites, terra::extract(elevation_PE, Per_Tahuayo_piloto_sites))



# get which are in and out
covs_Per_Tahuayo_sites$in_AP = st_intersects(covs_Per_Tahuayo_sites, AP_Tahuayo, sparse = FALSE)

covs_Per_Tahuayo_piloto_sites$in_AP = st_intersects(covs_Per_Tahuayo_piloto_sites, AP_Tahuayo, sparse = FALSE)




# make a map
mapview (elevation_PE, alpha=0.5) + 
  mapview (AP_Pacaya, color = "green", col.regions = "green", alpha = 0.5) +
  mapview (AP_Tahuayo, color = "green", col.regions = "green", alpha = 0.5) +
  mapview (covs_Per_Tahuayo_sites, zcol = "in_AP", col.regions =c("red","blue"), burst = TRUE) +   
  mapview (covs_Per_Tahuayo_piloto_sites, 
           zcol = "in_AP", 
           col.regions =c("blue"), 
           burst = TRUE) 

Yasuni data

Here we use the tables Ecu-13, Ecu-17, Ecu-18 y Ecu-20

Code
source("C:/CodigoR/WCS-CameraTrap/R/organiza_datos_v3.R")

AP_Yasuni <- read_sf("C:/CodigoR/Occu_APs/shp/Yasuni/WDPA_WDOECM_May2025_Public_186_shp-polygons.shp")

### Ecu 17, Ecu 18, ECU 20, Ecu 22  

# load data and make array_locID column
Ecu_13 <- loadproject("F:/WCS-CameraTrap/data/BDcorregidas/Ecuador/ECU-013.xlsx") |> mutate(array_locID=paste("Ecu_13", locationID, sep="_"))
28 cameras in Cameras. 
 28 cameras in Deployment. 
 28 deployments in Deployment. 
 28 points in Deployment. 
 28 cameras in Images. 
 28 points in Images. 
[1] "dates ok"
year: 2015 
 Jaguar_Design: no 
Code
Ecu_17 <- loadproject("F:/WCS-CameraTrap/data/BDcorregidas/Ecuador/ECU-017.xlsx") |> mutate(array_locID=paste("Ecu_17", locationID, sep="_"))
30 cameras in Cameras. 
 30 cameras in Deployment. 
 30 deployments in Deployment. 
 30 points in Deployment. 
 30 cameras in Images. 
 30 points in Images. 
[1] "dates ok"
year: 2015 
 Jaguar_Design: no 
Code
Ecu_18 <- loadproject("F:/WCS-CameraTrap/data/BDcorregidas/Ecuador/ECU-018.xlsx")|> mutate(array_locID=paste("Ecu_18", locationID, sep="_"))
30 cameras in Cameras. 
 30 cameras in Deployment. 
 30 deployments in Deployment. 
 30 points in Deployment. 
 30 cameras in Images. 
 30 points in Images. 
[1] "dates ok"
year: 2016 
 Jaguar_Design: no 
Code
Ecu_20 <- loadproject("F:/WCS-CameraTrap/data/BDcorregidas/Ecuador/ECU-020_Fix.xlsx")|> mutate(array_locID=paste("Ecu_20", locationID, sep="_"))
30 cameras in Cameras. 
 25 cameras in Deployment. 
 25 deployments in Deployment. 
 25 points in Deployment. 
 25 cameras in Images. 
 25 points in Images. 
[1] "dates ok"
year: 2018 
 Jaguar_Design: no 
Code
# Ecu_21 <- loadproject("F:/WCS-CameraTrap/data/BDcorregidas/Ecuador/ECU-021.xlsx") |> mutate(array_locID=paste("Ecu_14", locationID, sep="_"))
Ecu_22 <- loadproject("F:/WCS-CameraTrap/data/BDcorregidas/Ecuador/ECU-022.xlsx")|> mutate(array_locID=paste("Ecu_22", locationID, sep="_"))
30 cameras in Cameras. 
 30 cameras in Deployment. 
 30 deployments in Deployment. 
 30 points in Deployment. 
 30 cameras in Images. 
 30 points in Images. 
[1] "dates ok"
year: 2018 
 Jaguar_Design: no 
Code
# get sites
Ecu_13_sites <- get.sites("F:/WCS-CameraTrap/data/BDcorregidas/Ecuador/ECU-013.xlsx")
Ecu_17_sites <- get.sites("F:/WCS-CameraTrap/data/BDcorregidas/Ecuador/ECU-017.xlsx")
Ecu_18_sites <- get.sites("F:/WCS-CameraTrap/data/BDcorregidas/Ecuador/ECU-018.xlsx")
Ecu_20_sites <- get.sites("F:/WCS-CameraTrap/data/BDcorregidas/Ecuador/ECU-020_Fix.xlsx")
# Ecu_21_sites <- get.sites("F:/WCS-CameraTrap/data/BDcorregidas/Ecuador/ECU-021.xlsx")
Ecu_22_sites <- get.sites("F:/WCS-CameraTrap/data/BDcorregidas/Ecuador/ECU-022.xlsx")




# get elevation map
elevation_EC <- rast(get_elev_raster(Ecu_17_sites, z = 7)) #z =1-14
bb <-  st_as_sfc(st_bbox(elevation_EC)) # make bounding box 




# extract covs using points and add to _sites
covs_Ecu_13_sites <- cbind(Ecu_13_sites, terra::extract(elevation_EC, Ecu_13_sites))
covs_Ecu_17_sites <- cbind(Ecu_17_sites, terra::extract(elevation_EC, Ecu_17_sites))
covs_Ecu_18_sites <- cbind(Ecu_18_sites, terra::extract(elevation_EC, Ecu_18_sites))
covs_Ecu_20_sites <- cbind(Ecu_20_sites, terra::extract(elevation_EC, Ecu_20_sites))
#covs_Ecu_21_sites <- cbind(Ecu_21_sites, terra::extract(elevation_17, Ecu_21_sites))
covs_Ecu_22_sites <- cbind(Ecu_22_sites, terra::extract(elevation_EC, Ecu_22_sites))

# get which are in and out
covs_Ecu_13_sites$in_AP = st_intersects(covs_Ecu_13_sites, AP_Yasuni, sparse = FALSE)
covs_Ecu_17_sites$in_AP = st_intersects(covs_Ecu_17_sites, AP_Yasuni, sparse = FALSE)
covs_Ecu_18_sites$in_AP = st_intersects(covs_Ecu_18_sites, AP_Yasuni, sparse = FALSE)
covs_Ecu_20_sites$in_AP = st_intersects(covs_Ecu_20_sites, AP_Yasuni, sparse = FALSE)
#covs_Ecu_21_sites$in_AP = st_intersects(covs_Ecu_21_sites, AP_Yasuni, sparse = FALSE)
covs_Ecu_22_sites$in_AP = st_intersects(covs_Ecu_22_sites, AP_Yasuni, sparse = FALSE)

# covs_Ecu_16_sites$in_AP = st_intersects(covs_Ecu_16_sites, AP_Llanganates, sparse = FALSE)



# make a map
mapview (elevation_EC, alpha=0.5) + 
  mapview (AP_Yasuni, color = "green", col.regions = "green", alpha = 0.5) +
  mapview (covs_Ecu_13_sites, zcol = "in_AP", col.regions =c("red"), burst = TRUE) +
  mapview (covs_Ecu_17_sites, zcol = "in_AP", col.regions =c("red","blue"), burst = TRUE) +
  mapview (covs_Ecu_18_sites, zcol = "in_AP", col.regions =c("red"), burst = TRUE) +
  mapview (covs_Ecu_20_sites, zcol = "in_AP", col.regions =c("red","blue"), burst = TRUE) +
#  mapview (covs_Ecu_21_sites, zcol = "in_AP", col.regions =c("red"), burst = TRUE) +
  mapview (covs_Ecu_22_sites, zcol = "in_AP", burst = TRUE, col.regions = c("red") ) #+
Code
  # mapview (covs_Ecu_16_sites, zcol = "in_AP", burst = TRUE, col.regions =c("red","blue")) 

Camera trap operation data and detection history

Code
# Join 3 tables
# fix count in ECU 13, 17, 22,
Ecu_13$Count <- as.character(Ecu_13$Count)
Ecu_17$Count <- as.character(Ecu_17$Count)
Ecu_22$Count <- as.character(Ecu_22$Count)

Ecu_full <- Ecu_13 |> full_join(Ecu_17) |> 
                      full_join(Ecu_18) |> 
                      full_join(Ecu_20) |> 
#                      full_join(Ecu_21) |> 
                      full_join(Ecu_22)

###### Bolivia

Bol_full <- Bol_Madidi_1  |> 
   full_join(Bol_Madidi_2) |> 
   full_join(Bol_Madidi_3) |> 
   full_join(Bol_Madidi_4) 
# change camera Id by point. two cameras in one
Bol_full <- Bol_full |> mutate(Camera_Id=Point.x)


# Ecu_18$Count <- as.numeric(Ecu_18$Count)
Per_Tahuayo_piloto$Longitude <- as.numeric(Per_Tahuayo_piloto$Longitude)
Per_Tahuayo_piloto$Latitude <- as.numeric(Per_Tahuayo_piloto$Latitude)


Per_full <- Per_Tahuayo|> 
  full_join(Per_Tahuayo_piloto) #|> 
                      # full_join(Ecu_18) |> 
                      # full_join(Ecu_20) |> 
                      # full_join(Ecu_21) |> 
                      # full_join(Ecu_22)

# rename camera id
# Per_full$camid <- Per_full$`Camera_Id`

Ecu_full$Count <- as.numeric(Ecu_full$Count)
Per_full$Count <- as.numeric(Per_full$Count)

Ecu_full <- Ecu_full[,-47]
Col_18 <- Col_18[,-47]
##################################
data_full <- rbind(Ecu_full, 
                   Per_full,
                   Bol_full,
                   Col_18)


# fix date format
# 
# Formatting a Date object
data_full$start_date <- as.Date(data_full$"start_date", "%Y/%m/%d")
data_full$start_date <- format(data_full$start_date, "%Y-%m-%d")

data_full$end_date <- as.Date(data_full$"end_date", "%Y/%m/%d")
data_full$end_date <- format(data_full$end_date, "%Y-%m-%d")

data_full$eventDate <- as.Date(data_full$"Date_Time_Captured", "%Y/%m/%d")
data_full$eventDate <- format(data_full$eventDate, "%Y-%m-%d")

# Per_full$eventDateTime <- ymd_hms(paste(Per_full$"photo_date", Per_full$"photo_time", sep=" "))
data_full$eventDateTime <- ymd_hms(data_full$"Date_Time_Captured")

###############################
# remove duplicated cameras
################################
ind1 <- which(data_full$Camera_Id=="ECU-020-C0027")
ind2 <- which(data_full$Camera_Id=="ECU-020-C0006")
data_full <- data_full[-ind1,]
data_full <- data_full[-ind2,]


# filter 2021 and make uniques
CToperation  <- data_full |> dplyr::group_by(Camera_Id) |> #(array_locID) |> 
                           mutate(minStart=start_date, maxEnd=end_date) |> distinct(Longitude, Latitude, minStart, maxEnd) |> dplyr::ungroup()
# remove one duplicated
# View(CToperation)
# CToperation <- CToperation[-15,]

# M3-19: M3-21: M4-10: M4-12: M4-18:  From Madidi
CToperation <- CToperation[-c(379, 373, 365, 362, 357),]


CToperation[231,3] <- "-4.482831" #Latitude
CToperation[93,3] <- "-0.5548211" #Latitude
CToperation[93,2] <- "-76.48333"

### remove duplicated
# View(as.data.frame(table(CToperation$Latitude)))
# View(as.data.frame(table(CToperation$Longitude)))
# CToperation <- CToperation[-c(90,94),]
# Latitude -4.48283
# Latitude -0.554820969700813
# Longitude -76.4833290316164

# Generamos la matríz de operación de las cámaras

camop <- cameraOperation(CTtable= CToperation, # Tabla de operación
                         stationCol= "Camera_Id", # Columna que define la estación
                         setupCol= "minStart", #Columna fecha de colocación
                         retrievalCol= "maxEnd", #Columna fecha de retiro
                         #hasProblems= T, # Hubo fallos de cámaras
                         dateFormat= "%Y-%m-%d")#, #, # Formato de las fechas
                         #cameraCol="Camera_Id")
                         # sessionCol= "Year")

# Generar las historias de detección ---------------------------------------
## remove plroblem species

# Per_full$scientificName <- paste(Per_full$genus, Per_full$species, sep=" ")

#### remove NAs and setups
ind <- which(is.na(data_full$scientificName))
data_full <- data_full[-ind,]
# ind <- which(Per_full$scientificName=="NA NA")
# Per_full <- Per_full[-ind,]

# ind <- which(Per_full$scientificName=="Set up")
# Per_full <- Per_full[-ind,]
# 
# ind <- which(Per_full$scientificName=="Blank")
# Per_full <- Per_full[-ind,]
# 
# ind <- which(Per_full$scientificName=="Unidentifiable")
# Per_full <- Per_full[-ind,]


DetHist_list <- lapply(unique(data_full$scientificName), FUN = function(x) {
  detectionHistory(
    recordTable         = data_full, # abla de registros
    camOp                = camop, # Matriz de operación de cámaras
    stationCol           = "Camera_Id",
    speciesCol           = "scientificName",
    recordDateTimeCol    = "eventDateTime",
    recordDateTimeFormat  = "%Y-%m-%d %H:%M:%S",
    species              = x,     # la función reemplaza x por cada una de las especies
    occasionLength       = 7, # Colapso de las historias a 10 ías
    day1                 = "station", # "survey" a specific date, "station", #inicie en la fecha de cada survey
    datesAsOccasionNames = FALSE,
    includeEffort        = TRUE,
    scaleEffort          = FALSE,
    #unmarkedMultFrameInput=TRUE
    timeZone             = "America/Bogota" 
    )
  }
)

# names
names(DetHist_list) <- unique(data_full$scientificName)

# Finalmente creamos una lista nueva donde estén solo las historias de detección
ylist <- lapply(DetHist_list, FUN = function(x) x$detection_history)
# otra lista con effort scaled
efort <- lapply(DetHist_list, FUN = function(x) x$effort)

# number of observetions per sp, collapsed to 7 days
# lapply(ylist, sum, na.rm = TRUE)

Arrange spatial covariates

The standard EPSG code for a Lambert Azimuthal Equal-Area projection for South America is EPSG:10603 (WGS 84 / GLANCE South America). This projection is specifically designed for the continent and has its center located around the central meridian for the region. The units are meters.

FLII scores range from 0 (lowest integrity) to 10 (highest). Grantham discretized this range to define three broad illustrative categories: low (≤6.0); medium (>6.0 and <9.6); and high integrity (≥9.6).

Code
#transform coord data to Lambert Azimuthal Equal-Area
AP_Tahuayo_UTM <- st_transform(AP_Tahuayo, "EPSG:10603")
# Convert to LINESTRING
AP_Tahuayo_UTM_line <- st_cast(AP_Tahuayo_UTM,"MULTILINESTRING")# "LINESTRING")
#transform Yasuni to Lambert Azimuthal Equal-Area
AP_Yasuni_UTM <- st_transform(AP_Yasuni, "EPSG:10603")
# Convert to LINESTRING
AP_Yasuni_UTM_line <- st_cast(AP_Yasuni_UTM, "MULTILINESTRING")

#transform Yasuni to Lambert Azimuthal Equal-Area
AP_Madidi_UTM <- st_transform(Madidi_NP, "EPSG:10603")
# Convert to LINESTRING
AP_Madidi_UTM_line <- st_cast(AP_Madidi_UTM, "MULTILINESTRING")

AP_PlantasMed_UTM <- st_transform(AP_PlantasMed, "EPSG:10603")


############## sf AP Union 
AP_merged_sf_UTM <- st_union(AP_Tahuayo_UTM, 
                             AP_Yasuni_UTM,
                             AP_Madidi_UTM,
                             AP_PlantasMed_UTM)
Warning: attribute variables are assumed to be spatially constant throughout
all geometries
Code
AP_merged_line <- st_cast(AP_merged_sf_UTM, to="MULTILINESTRING")


# make sf() from data table
data_full_sf <- CToperation |> 
    st_as_sf(coords = c("Longitude", "Latitude"), 
              crs = 4326)

# gete elevation proint from AWS
elev_data <- elevatr::get_elev_point(locations = data_full_sf, src = "aws")
Mosaicing & Projecting
Note: Elevation units are in meters
Code
# extract elev and paste to table
data_full_sf$elev <- elev_data$elevation
str(data_full_sf$elev)
 num [1:428] 443 387 592 543 498 396 521 498 451 498 ...
Code
# extract in AP
data_full_sf$in_AP = as.factor(st_intersects(data_full_sf, AP_Tahuayo, sparse = FALSE))

in_AP <- as.numeric((st_drop_geometry(data_full_sf$in_AP)))

# extract FLII
data_full_sf$FLII <- terra::extract(FLII2016, data_full_sf)[,2]
str(data_full_sf$FLII)
 num [1:428] 9.57 9.65 9.69 9.66 9.92 ...
Code
# Replace all NAs with min flii in a numeric column
data_full_sf$FLII[is.na(data_full_sf$FLII)] <- min(data_full_sf$FLII, na.rm = TRUE)

# mapview(full_sites_14_15_16_sf, zcol = "in_AP", burst = TRUE)

# Transform coord to Lambert Azimuthal Equal-Area
data_full_sf_UTM <- st_transform(data_full_sf, "EPSG:10603")


coords <- st_coordinates(data_full_sf_UTM)
#str(coords)

#### fix duplicated coord
# -1840583.53296873  en x
# -1842515.36969736  en x
# 1547741.44311964  en y
# 1541202.24796904 en y

# which(coords[,1]=="-1840583.53296873")
# which(coords[,1]=="-1842515.36969736")
# which(coords[,2]=="1547741.44311964")
# which(coords[,2]=="1541202.24796904")

# make Ecu_14_15_16 an sf object
#    cam_sf <- st_as_sf(Ecu_14_15_16, coords = c("lon","lat"))   #crs="EPSG:4326")
    #--- set CRS ---#
#    st_crs(cam_sf) <- 4326




# Calculate the distance
#multiplic <- full_sites_14_15_16_sf_UTM |> mutate(multiplic= as.numeric(in_AP)) 
multiplic=ifelse(data_full_sf_UTM$in_AP=="TRUE",-1,1)
data_full_sf_UTM$border_dist <- as.numeric(st_distance(data_full_sf_UTM, AP_Tahuayo_UTM_line) * multiplic )
# print(border_dist)

# convert true false to inside 1, outside 0
data_full_sf_UTM <- data_full_sf_UTM |>
  mutate(in_AP = case_when(
    str_detect(in_AP, "TRUE") ~ 1, # "inside_AP",
    str_detect(in_AP, "FALSE") ~ 0 #"outside_AP"
  )) |> mutate(in_AP=as.factor(in_AP))


hist(data_full_sf_UTM$border_dist)

Prepare the model

TipData in a 3D array

The data must be placed in a three-dimensional array with dimensions corresponding to species, sites, and replicates in that order.

The function sfMsPGOcc

Fits multi-species spatial occupancy models with species correlations (i.e., a spatially-explicit joint species distribution model with imperfect detection). We use Polya-Gamma latent variables and a spatial factor modeling approach. Currently, models are implemented using a Nearest Neighbor Gaussian Process.

Code
# Detection-nondetection data ---------
# Species of interest, can select individually
# curr.sp <- sort(unique(Ecu_14_15_16$.id))# c('BAWW', 'BLJA', 'GCFL')
# sort(names(DetHist_list))
selected.sp <-  c(
"Atelocynus microtis" ,
"Coendou prehensilis" ,
"Cuniculus paca",           
"Dasyprocta fuliginosa",      
"Dasypus sp." ,    
"Didelphis marsupialis",    
"Eira barbara",  
"Herpailurus yagouaroundi",
"Leopardus pardalis",    
"Leopardus tigrinus" ,
"Leopardus wiedii",         
"Mazama americana",
"Mazama nemorivaga",
#"Mitu tuberosum" ,
"Myoprocta pratti",
"Myrmecophaga tridactyla",
"Nasua nasua" ,
# "Mazama americana",         
# "Myotis myotis",           
# "Nasua narica",             
# "Odocoileus virginianus",   
"Panthera onca" ,
# "Procyon cancrivorus" ,
"Pecari tajacu",    
#"Penelope jacquacu" ,
"Priodontes maximus" ,
"Procyon cancrivorous",      
# "Psophia leucoptera",
"Puma concolor" ,
"Puma yagouaroundi",        
# "Rattus rattus" ,
# "Roedor sp.",
# "Sciurus sp.",       
# "Sus scrofa",               
# "Sylvilagus brasiliensis",  
"Tamandua tetradactyla",   
"Tapirus terrestris",
"Tayassu pecari"
#"Tinamus major"            
              )

# y.msom <- y[which(sp.codes %in% selected.sp), , ]
# str(y.msom)

# Use selection
y.selected <- ylist[selected.sp]   

#### three-dimensional array with dimensions corresponding to species, sites, and replicates

# 1. Load the abind library to make arrays easily 
library(abind)
my_array_abind <- abind(y.selected, # start from list
                        along = 3, # 3D array
                        use.first.dimnames=TRUE) # keep names

# Transpose the array to have:
# species, sites, and sampling occasions in that order
# The new order is (3rd dim, 1st dim, 2nd dim)
transposed_array <- aperm(my_array_abind, c(3, 1, 2))

#### site covs
sitecovs <- as.data.frame(st_drop_geometry(
                    data_full_sf_UTM[,5:8]))

 sitecovs[, 1] <- as.vector((sitecovs[,1]))   # scale numeric covariates
 sitecovs[, 1] <- as.numeric((sitecovs[,2]))   # scale numeric covariates
 sitecovs[, 3] <- as.vector((sitecovs[,3]))   # scale numeric covariates
 sitecovs[, 4] <- as.vector((sitecovs[,4]))   # scale numeric covariates
 # sitecovs$fact <- factor(c("A", "A", "B"))    # categorical covariate

names(sitecovs) <- c("elev", "in_AP", "FLII", "border_dist")

# check consistancy equal number of spatial covariates and rows in data
# identical(nrow(ylist[[1]]), nrow(covars)) 

# Base de datos para los análisis -----------------------------------------

# match the names to "y"  "occ.covs" "det.covs" "coords" 
data_list <- list(y = transposed_array, # Historias de detección
                  occ.covs = sitecovs, # covs de sitio
                  det.covs  = list(effort = DetHist_list[[1]]$effort),
                  coords = st_coordinates(data_full_sf_UTM)
                  )  # agregamos el esfuerzo de muestreo como covariable de observación

Running the model

We let spOccupancy set the initial values by default based on the prior distributions.

Code
# Running the model

# 3. 1 Modelo multi-especie  -----------------------------------------

# 2. Model fitting --------------------------------------------------------
# Fit a non-spatial, multi-species occupancy model. 
out <- msPGOcc(occ.formula = ~ scale(elev) +
                               scale(border_dist) + 
                               scale(FLII) , 
               det.formula = ~ scale(effort) , # Ordinal.day + I(Ordinal.day^2) + Year
                 data = data_list, 
                 n.samples = 6000, 
                 n.thin = 10, 
                 n.burn = 1000, 
                 n.chains = 3,
                 n.report = 1000);beep(sound = 4)
----------------------------------------
    Preparing to run the model
----------------------------------------
No prior specified for beta.comm.normal.
Setting prior mean to 0 and prior variance to 2.72
No prior specified for alpha.comm.normal.
Setting prior mean to 0 and prior variance to 2.72
No prior specified for tau.sq.beta.ig.
Setting prior shape to 0.1 and prior scale to 0.1
No prior specified for tau.sq.alpha.ig.
Setting prior shape to 0.1 and prior scale to 0.1
z is not specified in initial values.
Setting initial values based on observed data
beta.comm is not specified in initial values.
Setting initial values to random values from the prior distribution
alpha.comm is not specified in initial values.
Setting initial values to random values from the prior distribution
tau.sq.beta is not specified in initial values.
Setting initial values to random values between 0.5 and 10
tau.sq.alpha is not specified in initial values.
Setting to initial values to random values between 0.5 and 10
beta is not specified in initial values.
Setting initial values to random values from the community-level normal distribution
alpha is not specified in initial values.
Setting initial values to random values from the community-level normal distribution
----------------------------------------
    Model description
----------------------------------------
Multi-species Occupancy Model with Polya-Gamma latent
variable fit with 428 sites and 25 species.

Samples per Chain: 6000 
Burn-in: 1000 
Thinning Rate: 10 
Number of Chains: 3 
Total Posterior Samples: 1500 

Source compiled with OpenMP support and model fit using 1 thread(s).

----------------------------------------
    Chain 1
----------------------------------------
Sampling ... 
Sampled: 1000 of 6000, 16.67%
-------------------------------------------------
Sampled: 2000 of 6000, 33.33%
-------------------------------------------------
Sampled: 3000 of 6000, 50.00%
-------------------------------------------------
Sampled: 4000 of 6000, 66.67%
-------------------------------------------------
Sampled: 5000 of 6000, 83.33%
-------------------------------------------------
Sampled: 6000 of 6000, 100.00%
----------------------------------------
    Chain 2
----------------------------------------
Sampling ... 
Sampled: 1000 of 6000, 16.67%
-------------------------------------------------
Sampled: 2000 of 6000, 33.33%
-------------------------------------------------
Sampled: 3000 of 6000, 50.00%
-------------------------------------------------
Sampled: 4000 of 6000, 66.67%
-------------------------------------------------
Sampled: 5000 of 6000, 83.33%
-------------------------------------------------
Sampled: 6000 of 6000, 100.00%
----------------------------------------
    Chain 3
----------------------------------------
Sampling ... 
Sampled: 1000 of 6000, 16.67%
-------------------------------------------------
Sampled: 2000 of 6000, 33.33%
-------------------------------------------------
Sampled: 3000 of 6000, 50.00%
-------------------------------------------------
Sampled: 4000 of 6000, 66.67%
-------------------------------------------------
Sampled: 5000 of 6000, 83.33%
-------------------------------------------------
Sampled: 6000 of 6000, 100.00%
Code
# Fit a non-spatial, Latent Factor Multi-Species Occupancy Model. 
  out.lfMs <- msPGOcc(occ.formula = ~ scale(elev) +
                               scale(border_dist) + 
                               scale(FLII) , 
               det.formula = ~ scale(effort), # Ordinal.day + I(Ordinal.day^2) + Year
                 data = data_list, 
                 n.omp.threads = 6,
                 n.samples = 6000, 
                 n.factors = 5, # balance of rare sp. and run time
                 n.thin = 10, 
                 n.burn = 1000, 
                 n.chains = 3,
                 n.report = 1000);beep(sound = 4)
----------------------------------------
    Preparing to run the model
----------------------------------------
Warning in msPGOcc(occ.formula = ~scale(elev) + scale(border_dist) +
scale(FLII), : 'n.factors' is not an argument
No prior specified for beta.comm.normal.
Setting prior mean to 0 and prior variance to 2.72
No prior specified for alpha.comm.normal.
Setting prior mean to 0 and prior variance to 2.72
No prior specified for tau.sq.beta.ig.
Setting prior shape to 0.1 and prior scale to 0.1
No prior specified for tau.sq.alpha.ig.
Setting prior shape to 0.1 and prior scale to 0.1
z is not specified in initial values.
Setting initial values based on observed data
beta.comm is not specified in initial values.
Setting initial values to random values from the prior distribution
alpha.comm is not specified in initial values.
Setting initial values to random values from the prior distribution
tau.sq.beta is not specified in initial values.
Setting initial values to random values between 0.5 and 10
tau.sq.alpha is not specified in initial values.
Setting to initial values to random values between 0.5 and 10
beta is not specified in initial values.
Setting initial values to random values from the community-level normal distribution
alpha is not specified in initial values.
Setting initial values to random values from the community-level normal distribution
----------------------------------------
    Model description
----------------------------------------
Multi-species Occupancy Model with Polya-Gamma latent
variable fit with 428 sites and 25 species.

Samples per Chain: 6000 
Burn-in: 1000 
Thinning Rate: 10 
Number of Chains: 3 
Total Posterior Samples: 1500 

Source compiled with OpenMP support and model fit using 6 thread(s).

----------------------------------------
    Chain 1
----------------------------------------
Sampling ... 
Sampled: 1000 of 6000, 16.67%
-------------------------------------------------
Sampled: 2000 of 6000, 33.33%
-------------------------------------------------
Sampled: 3000 of 6000, 50.00%
-------------------------------------------------
Sampled: 4000 of 6000, 66.67%
-------------------------------------------------
Sampled: 5000 of 6000, 83.33%
-------------------------------------------------
Sampled: 6000 of 6000, 100.00%
----------------------------------------
    Chain 2
----------------------------------------
Sampling ... 
Sampled: 1000 of 6000, 16.67%
-------------------------------------------------
Sampled: 2000 of 6000, 33.33%
-------------------------------------------------
Sampled: 3000 of 6000, 50.00%
-------------------------------------------------
Sampled: 4000 of 6000, 66.67%
-------------------------------------------------
Sampled: 5000 of 6000, 83.33%
-------------------------------------------------
Sampled: 6000 of 6000, 100.00%
----------------------------------------
    Chain 3
----------------------------------------
Sampling ... 
Sampled: 1000 of 6000, 16.67%
-------------------------------------------------
Sampled: 2000 of 6000, 33.33%
-------------------------------------------------
Sampled: 3000 of 6000, 50.00%
-------------------------------------------------
Sampled: 4000 of 6000, 66.67%
-------------------------------------------------
Sampled: 5000 of 6000, 83.33%
-------------------------------------------------
Sampled: 6000 of 6000, 100.00%
Code
summary(out, level = 'community')

Call:
msPGOcc(occ.formula = ~scale(elev) + scale(border_dist) + scale(FLII), 
    det.formula = ~scale(effort), data = data_list, n.samples = 6000, 
    n.report = 1000, n.burn = 1000, n.thin = 10, n.chains = 3)

Samples per Chain: 6000
Burn-in: 1000
Thinning Rate: 10
Number of Chains: 3
Total Posterior Samples: 1500
Run Time (min): 7.194

----------------------------------------
    Community Level
----------------------------------------
Occurrence Means (logit scale): 
                      Mean     SD    2.5%     50%   97.5%   Rhat  ESS
(Intercept)        -1.9989 0.3610 -2.7013 -2.0081 -1.2564 1.0327  394
scale(elev)        -0.1799 0.1238 -0.4193 -0.1780  0.0734 0.9994  934
scale(border_dist) -0.3663 0.2178 -0.7892 -0.3635  0.0553 1.0005 1047
scale(FLII)         0.4476 0.1907  0.0894  0.4425  0.8499 1.0038  745

Occurrence Variances (logit scale): 
                     Mean     SD   2.5%    50%  97.5%   Rhat ESS
(Intercept)        2.7795 1.0549 1.3005 2.5494 5.3226 1.0588 161
scale(elev)        0.2559 0.1191 0.1071 0.2294 0.5436 1.0145 991
scale(border_dist) 1.0008 0.4134 0.4484 0.9310 2.0193 1.0306 537
scale(FLII)        0.5658 0.2902 0.2031 0.5049 1.2841 1.0161 443

Detection Means (logit scale): 
                 Mean     SD    2.5%     50%   97.5%   Rhat  ESS
(Intercept)   -2.2353 0.2703 -2.8274 -2.2257 -1.7444 1.0426  502
scale(effort)  0.3975 0.0667  0.2619  0.3944  0.5320 0.9995 1178

Detection Variances (logit scale): 
                Mean     SD   2.5%    50%  97.5%   Rhat  ESS
(Intercept)   1.4161 0.6351 0.6087 1.2964 2.9403 1.0614  139
scale(effort) 0.0480 0.0278 0.0175 0.0410 0.1191 1.0062 1195
Code
# Fit a Spatial Factor Multi-Species Occupancy Model
# latent spatial factors.
tictoc::tic()
  out.sp <- sfMsPGOcc(occ.formula = ~ scale(elev) +
                               scale(border_dist) + 
                               scale(FLII) , 
               det.formula = ~ scale(effort), # Ordinal.day + I(Ordinal.day^2) + Year,            
                      data = data_list, 
                      n.omp.threads = 6,
                      n.batch = 600, 
                      batch.length = 25, # iter=600*25
                      n.thin = 10, 
                      n.burn = 5000, 
                      n.chains = 3,
                      NNGP = TRUE,
                      n.factors = 5, # balance of rare sp. and run time
                      n.neighbors = 15,
                      cov.model = 'exponential',
                      n.report = 1000);beep(sound = 4)
----------------------------------------
    Preparing to run the model
----------------------------------------
No prior specified for beta.comm.normal.
Setting prior mean to 0 and prior variance to 2.72
No prior specified for alpha.comm.normal.
Setting prior mean to 0 and prior variance to 2.72
No prior specified for tau.sq.beta.
Using an inverse-Gamma prior with prior shape 0.1 and prior scale 0.1
No prior specified for tau.sq.alpha.
Using an inverse-Gamma prior with prior shape 0.1 and prior scale 0.1
No prior specified for phi.unif.
Setting uniform bounds based on the range of observed spatial coordinates.
z is not specified in initial values.
Setting initial values based on observed data
beta.comm is not specified in initial values.
Setting initial values to random values from the prior distribution
alpha.comm is not specified in initial values.
Setting initial values to random values from the prior distribution
tau.sq.beta is not specified in initial values.
Setting initial values to random values between 0.5 and 10
tau.sq.alpha is not specified in initial values.
Setting initial values to random values between 0.5 and 10
beta is not specified in initial values.
Setting initial values to random values from the community-level normal distribution
alpha is not specified in initial values.
Setting initial values to random values from the community-level normal distribution
phi is not specified in initial values.
Setting initial value to random values from the prior distribution
lambda is not specified in initial values.
Setting initial values of the lower triangle to 0
----------------------------------------
    Building the neighbor list
----------------------------------------
----------------------------------------
Building the neighbors of neighbors list
----------------------------------------
----------------------------------------
    Model description
----------------------------------------
Spatial Factor NNGP Multi-species Occupancy Model with Polya-Gamma latent
variable fit with 428 sites and 25 species.

Samples per chain: 15000 (600 batches of length 25)
Burn-in: 5000 
Thinning Rate: 10 
Number of Chains: 3 
Total Posterior Samples: 3000 

Using the exponential spatial correlation model.

Using 5 latent spatial factors.
Using 15 nearest neighbors.

Source compiled with OpenMP support and model fit using 6 thread(s).

Adaptive Metropolis with target acceptance rate: 43.0
----------------------------------------
    Chain 1
----------------------------------------
Sampling ... 
Batch: 600 of 600, 100.00%
----------------------------------------
    Chain 2
----------------------------------------
Sampling ... 
Batch: 600 of 600, 100.00%
----------------------------------------
    Chain 3
----------------------------------------
Sampling ... 
Batch: 600 of 600, 100.00%
Code
tictoc::toc()
1800.95 sec elapsed
Code
#########################
tictoc::tic()
  out.sp.fac <- sfMsPGOcc(occ.formula = ~ scale(elev) +
                               factor(in_AP) + 
                               scale(FLII) , 
               det.formula = ~ scale(effort), # Ordinal.day + I(Ordinal.day^2) + Year,            
                      data = data_list, 
                      n.omp.threads = 6,
                      n.batch = 600, 
                      batch.length = 25, # iter=600*25
                      n.thin = 10, 
                      n.burn = 5000, 
                      n.chains = 3,
                      NNGP = TRUE,
                      n.factors = 5, # balance of rare sp. and run time
                      n.neighbors = 15,
                      cov.model = 'exponential',
                      n.report = 1000);beep(sound = 4)
----------------------------------------
    Preparing to run the model
----------------------------------------
No prior specified for beta.comm.normal.
Setting prior mean to 0 and prior variance to 2.72
No prior specified for alpha.comm.normal.
Setting prior mean to 0 and prior variance to 2.72
No prior specified for tau.sq.beta.
Using an inverse-Gamma prior with prior shape 0.1 and prior scale 0.1
No prior specified for tau.sq.alpha.
Using an inverse-Gamma prior with prior shape 0.1 and prior scale 0.1
No prior specified for phi.unif.
Setting uniform bounds based on the range of observed spatial coordinates.
z is not specified in initial values.
Setting initial values based on observed data
beta.comm is not specified in initial values.
Setting initial values to random values from the prior distribution
alpha.comm is not specified in initial values.
Setting initial values to random values from the prior distribution
tau.sq.beta is not specified in initial values.
Setting initial values to random values between 0.5 and 10
tau.sq.alpha is not specified in initial values.
Setting initial values to random values between 0.5 and 10
beta is not specified in initial values.
Setting initial values to random values from the community-level normal distribution
alpha is not specified in initial values.
Setting initial values to random values from the community-level normal distribution
phi is not specified in initial values.
Setting initial value to random values from the prior distribution
lambda is not specified in initial values.
Setting initial values of the lower triangle to 0
----------------------------------------
    Building the neighbor list
----------------------------------------
----------------------------------------
Building the neighbors of neighbors list
----------------------------------------
----------------------------------------
    Model description
----------------------------------------
Spatial Factor NNGP Multi-species Occupancy Model with Polya-Gamma latent
variable fit with 428 sites and 25 species.

Samples per chain: 15000 (600 batches of length 25)
Burn-in: 5000 
Thinning Rate: 10 
Number of Chains: 3 
Total Posterior Samples: 3000 

Using the exponential spatial correlation model.

Using 5 latent spatial factors.
Using 15 nearest neighbors.

Source compiled with OpenMP support and model fit using 6 thread(s).

Adaptive Metropolis with target acceptance rate: 43.0
----------------------------------------
    Chain 1
----------------------------------------
Sampling ... 
Batch: 600 of 600, 100.00%
----------------------------------------
    Chain 2
----------------------------------------
Sampling ... 
Batch: 600 of 600, 100.00%
----------------------------------------
    Chain 3
----------------------------------------
Sampling ... 
Batch: 600 of 600, 100.00%
Code
tictoc::toc()
1795.69 sec elapsed
Code
# tictoc::tic()
#   out.sp.gaus <- sfMsPGOcc(occ.formula = ~ scale(border_dist) , 
#                            det.formula = ~ scale(effort), # Ordinal.day + I(Ordinal.day^2) + Year,              
#                            data = data_list, 
#                            n.batch = 400, 
#                            batch.length = 25,
#                            n.thin = 5, 
#                            n.burn = 5000, 
#                            n.chains = 1,
#                            NNGP = TRUE,
#                            n.factors = 5,
#                            n.neighbors = 15,
#                            cov.model = 'gaussian',
#                            n.report = 100);beep(sound = 4)
# tictoc::toc()

# save the results to not run again
# save(out, file="C:/CodigoR/Occu_APs_all/blog/2025-10-15-analysis/result/result_2.R") # guardamos los resultados para no correr de nuevo
# save the results to not run again
# save(out.sp, file="C:/CodigoR/Occu_APs_all/blog/2025-10-15-analysis/result/sp_result_2.R") # guardamos los resultados para no correr de nuevo

# save(out.lfMs, file="C:/CodigoR/Occu_APs_all/blog/2025-10-15-analysis/result/lfms_result_2.R") # guardamos los resultados para no correr de nuevo


# load("C:/CodigoR/Occu_APs_all/blog/2025-10-15-analysis/result/sp_result_2.R")
# summary(fit.commu)

Model validation

We next perform a posterior predictive check using the Freeman-Tukey statistic grouping the data by sites. We summarize the posterior predictive check with the summary() function, which reports a Bayesian p-value. A Bayesian p-value that hovers around 0.5 indicates adequate model fit, while values less than 0.1 or greater than 0.9 suggest our model does not fit the data well (Hobbs and Hooten 2015). As always with a simulation-based analysis using MCMC, you will get numerically slightly different values.

Code
# 3. Model validation -----------------------------------------------------
# Perform a posterior predictive check to assess model fit. 
ppc.out <- ppcOcc(out, fit.stat = 'freeman-tukey', 
                  group = 1)
ppc.out.lfMs <- ppcOcc(out.lfMs, fit.stat = 'freeman-tukey', 
                  group = 1)
ppc.out.sp <- ppcOcc(out.sp, fit.stat = 'freeman-tukey',
                     group = 1)

# Calculate a Bayesian p-value as a simple measure of Goodness of Fit.
# Bayesian p-values between 0.1 and 0.9 indicate adequate model fit. 
summary(ppc.out)

Call:
ppcOcc(object = out, fit.stat = "freeman-tukey", group = 1)

Samples per Chain: 6000
Burn-in: 1000
Thinning Rate: 10
Number of Chains: 3
Total Posterior Samples: 1500

----------------------------------------
    Community Level
----------------------------------------
Bayesian p-value:  0.3685 

----------------------------------------
    Species Level
----------------------------------------
Atelocynus microtis Bayesian p-value: 0.3367
Coendou prehensilis Bayesian p-value: 0.6347
Cuniculus paca Bayesian p-value: 0
Dasyprocta fuliginosa Bayesian p-value: 0.0013
Dasypus sp. Bayesian p-value: 0.1493
Didelphis marsupialis Bayesian p-value: 0.4647
Eira barbara Bayesian p-value: 0.9207
Herpailurus yagouaroundi Bayesian p-value: 0.5413
Leopardus pardalis Bayesian p-value: 0.3447
Leopardus tigrinus Bayesian p-value: 0.2013
Leopardus wiedii Bayesian p-value: 0.364
Mazama americana Bayesian p-value: 0.058
Mazama nemorivaga Bayesian p-value: 0.5067
Myoprocta pratti Bayesian p-value: 0.1467
Myrmecophaga tridactyla Bayesian p-value: 0.3953
Nasua nasua Bayesian p-value: 0.5153
Panthera onca Bayesian p-value: 0.3847
Pecari tajacu Bayesian p-value: 0.1713
Priodontes maximus Bayesian p-value: 0.516
Procyon cancrivorous Bayesian p-value: 0.5193
Puma concolor Bayesian p-value: 0.386
Puma yagouaroundi Bayesian p-value: 0.678
Tamandua tetradactyla Bayesian p-value: 0.676
Tapirus terrestris Bayesian p-value: 0.0693
Tayassu pecari Bayesian p-value: 0.232
Fit statistic:  freeman-tukey 
Code
summary(ppc.out.lfMs)

Call:
ppcOcc(object = out.lfMs, fit.stat = "freeman-tukey", group = 1)

Samples per Chain: 6000
Burn-in: 1000
Thinning Rate: 10
Number of Chains: 3
Total Posterior Samples: 1500

----------------------------------------
    Community Level
----------------------------------------
Bayesian p-value:  0.3661 

----------------------------------------
    Species Level
----------------------------------------
Atelocynus microtis Bayesian p-value: 0.312
Coendou prehensilis Bayesian p-value: 0.66
Cuniculus paca Bayesian p-value: 0
Dasyprocta fuliginosa Bayesian p-value: 0.002
Dasypus sp. Bayesian p-value: 0.1653
Didelphis marsupialis Bayesian p-value: 0.4587
Eira barbara Bayesian p-value: 0.9127
Herpailurus yagouaroundi Bayesian p-value: 0.53
Leopardus pardalis Bayesian p-value: 0.336
Leopardus tigrinus Bayesian p-value: 0.2227
Leopardus wiedii Bayesian p-value: 0.3587
Mazama americana Bayesian p-value: 0.08
Mazama nemorivaga Bayesian p-value: 0.4813
Myoprocta pratti Bayesian p-value: 0.1393
Myrmecophaga tridactyla Bayesian p-value: 0.3587
Nasua nasua Bayesian p-value: 0.5507
Panthera onca Bayesian p-value: 0.3833
Pecari tajacu Bayesian p-value: 0.1727
Priodontes maximus Bayesian p-value: 0.498
Procyon cancrivorous Bayesian p-value: 0.5233
Puma concolor Bayesian p-value: 0.3867
Puma yagouaroundi Bayesian p-value: 0.6753
Tamandua tetradactyla Bayesian p-value: 0.644
Tapirus terrestris Bayesian p-value: 0.078
Tayassu pecari Bayesian p-value: 0.224
Fit statistic:  freeman-tukey 
Code
summary(ppc.out.sp)

Call:
ppcOcc(object = out.sp, fit.stat = "freeman-tukey", group = 1)

Samples per Chain: 15000
Burn-in: 5000
Thinning Rate: 10
Number of Chains: 3
Total Posterior Samples: 3000

----------------------------------------
    Community Level
----------------------------------------
Bayesian p-value:  0.3537 

----------------------------------------
    Species Level
----------------------------------------
Atelocynus microtis Bayesian p-value: 0.33
Coendou prehensilis Bayesian p-value: 0.631
Cuniculus paca Bayesian p-value: 0
Dasyprocta fuliginosa Bayesian p-value: 0.0027
Dasypus sp. Bayesian p-value: 0.158
Didelphis marsupialis Bayesian p-value: 0.418
Eira barbara Bayesian p-value: 0.9497
Herpailurus yagouaroundi Bayesian p-value: 0.569
Leopardus pardalis Bayesian p-value: 0.2787
Leopardus tigrinus Bayesian p-value: 0.2447
Leopardus wiedii Bayesian p-value: 0.3077
Mazama americana Bayesian p-value: 0.0463
Mazama nemorivaga Bayesian p-value: 0.6223
Myoprocta pratti Bayesian p-value: 0.1463
Myrmecophaga tridactyla Bayesian p-value: 0.1767
Nasua nasua Bayesian p-value: 0.4403
Panthera onca Bayesian p-value: 0.4193
Pecari tajacu Bayesian p-value: 0.1733
Priodontes maximus Bayesian p-value: 0.5257
Procyon cancrivorous Bayesian p-value: 0.4837
Puma concolor Bayesian p-value: 0.351
Puma yagouaroundi Bayesian p-value: 0.688
Tamandua tetradactyla Bayesian p-value: 0.62
Tapirus terrestris Bayesian p-value: 0.0527
Tayassu pecari Bayesian p-value: 0.2087
Fit statistic:  freeman-tukey 
Important

A Bayesian p-value that around 0.5 indicates adequate model fit, while values less than 0.1 or greater than 0.9 suggest our model does not fit the data well.

Fit plot

Code
#### fit plot
ppc.df <- data.frame(fit = ppc.out.sp$fit.y, 
                     fit.rep = ppc.out.sp$fit.y.rep, 
                     color = 'lightskyblue1')

ppc.df$color[ppc.df$fit.rep.1 > ppc.df$fit.1] <- 'lightsalmon'
plot(ppc.df$fit.1, ppc.df$fit.rep.1, bg = ppc.df$color, pch = 21, 
     ylab = 'Fit', xlab = 'True')
lines(ppc.df$fit.1, ppc.df$fit.1, col = 'black')

The most symmetrical, better fit!

The most symmetrical, better fit!

Model comparison

Code
# 4. Model comparison -----------------------------------------------------
# Compute Widely Applicable Information Criterion (WAIC)
# Lower values indicate better model fit. 
waicOcc(out)
     elpd        pD      WAIC 
-8225.109   118.086 16686.389 
Code
waicOcc(out.lfMs)
      elpd         pD       WAIC 
-8224.7917   118.5652 16686.7138 
Code
waicOcc(out.sp)
      elpd         pD       WAIC 
-7280.6104   739.5953 16040.4113 
Code
waicOcc(out.sp.fac)
      elpd         pD       WAIC 
-7333.4507   712.2107 16091.3228 

Model comparison as table

Code
# Here we summarize the spatial factor loadings
# summary(out.sp$lambda.samples)

# Resultados --------------------------------------------------------------
# Extraemos lo tabla de valores estimados
modresult <- cbind(as.data.frame(waicOcc(out)),
                  as.data.frame(waicOcc(out.sp)),
                  as.data.frame(waicOcc(out.lfMs)),
                  as.data.frame(waicOcc(out.sp.fac))
                  #as.data.frame(waicOcc(out.sp.gaus))
                  )
# View(modresult)
modresult_sorted <- as.data.frame(t(modresult)) |> 
  arrange(WAIC) # sort by
  
DT::datatable(modresult_sorted)
Important

Lower values in WAIC indicate better model fit.

WAIC is the Widely Applicable Information Criterion (Watanabe 2010).

Posterior Summary

Code
# 5. Posterior summaries --------------------------------------------------
# Concise summary of main parameter estimates
summary(out.sp, level = 'community')

Call:
sfMsPGOcc(occ.formula = ~scale(elev) + scale(border_dist) + scale(FLII), 
    det.formula = ~scale(effort), data = data_list, cov.model = "exponential", 
    NNGP = TRUE, n.neighbors = 15, n.factors = 5, n.batch = 600, 
    batch.length = 25, n.omp.threads = 6, n.report = 1000, n.burn = 5000, 
    n.thin = 10, n.chains = 3)

Samples per Chain: 15000
Burn-in: 5000
Thinning Rate: 10
Number of Chains: 3
Total Posterior Samples: 3000
Run Time (min): 30.0153

----------------------------------------
    Community Level
----------------------------------------
Occurrence Means (logit scale): 
                      Mean     SD    2.5%     50%   97.5%   Rhat ESS
(Intercept)        -2.7611 0.4954 -3.7565 -2.7599 -1.7782 1.0807 277
scale(elev)        -0.0372 0.1944 -0.4196 -0.0382  0.3414 1.1540 850
scale(border_dist) -0.4589 0.3319 -1.1346 -0.4491  0.1698 1.0799 895
scale(FLII)         0.5170 0.2460  0.0567  0.5049  1.0430 1.0122 885

Occurrence Variances (logit scale): 
                     Mean     SD   2.5%    50%   97.5%   Rhat  ESS
(Intercept)        5.1766 2.1408 2.2988 4.7363 10.7314 1.1570  346
scale(elev)        0.4771 0.2324 0.1718 0.4252  1.0559 1.0385 1145
scale(border_dist) 2.0565 0.8866 0.8445 1.8913  4.2675 1.0374  438
scale(FLII)        0.9112 0.4685 0.3170 0.8046  2.0406 1.0008  496

Detection Means (logit scale): 
                 Mean     SD    2.5%     50%   97.5%   Rhat  ESS
(Intercept)   -2.2605 0.2700 -2.8176 -2.2461 -1.7495 1.1593  329
scale(effort)  0.3955 0.0653  0.2704  0.3935  0.5316 1.0033 2257

Detection Variances (logit scale): 
                Mean     SD   2.5%    50%  97.5%   Rhat  ESS
(Intercept)   1.5039 0.6865 0.6342 1.3541 3.2076 1.2673  106
scale(effort) 0.0454 0.0261 0.0160 0.0392 0.1110 1.0008 2566

----------------------------------------
    Spatial Covariance
----------------------------------------
         Mean     SD   2.5%     50%   97.5%   Rhat  ESS
phi-1  8.9235 9.0901 0.0000  6.4319 26.4508 4.0224    8
phi-2 13.2789 7.9866 0.4498 13.0502 26.7775 1.0019 3000
phi-3  4.7494 8.0864 0.0000  0.0000 26.0644 9.1416    3
phi-4 13.5261 7.9445 0.6362 13.4657 26.8778 1.0001 2677
phi-5 13.6591 8.0601 0.6792 13.3846 26.8695 1.0007 3000
Code
summary(out.sp, level = 'species')

Call:
sfMsPGOcc(occ.formula = ~scale(elev) + scale(border_dist) + scale(FLII), 
    det.formula = ~scale(effort), data = data_list, cov.model = "exponential", 
    NNGP = TRUE, n.neighbors = 15, n.factors = 5, n.batch = 600, 
    batch.length = 25, n.omp.threads = 6, n.report = 1000, n.burn = 5000, 
    n.thin = 10, n.chains = 3)

Samples per Chain: 15000
Burn-in: 5000
Thinning Rate: 10
Number of Chains: 3
Total Posterior Samples: 3000
Run Time (min): 30.0153

----------------------------------------
    Species Level
----------------------------------------
Occurrence (logit scale): 
                                               Mean     SD    2.5%     50%
(Intercept)-Atelocynus microtis             -2.5173 1.0719 -3.9337 -2.7944
(Intercept)-Coendou prehensilis             -5.4512 1.7234 -8.7614 -5.5698
(Intercept)-Cuniculus paca                   0.1873 0.4071 -0.5579  0.1865
(Intercept)-Dasyprocta fuliginosa           -1.3041 1.1359 -3.5264 -1.2099
(Intercept)-Dasypus sp.                     -5.2100 0.9039 -7.1616 -5.1349
(Intercept)-Didelphis marsupialis           -3.9183 0.8388 -5.6388 -3.8827
(Intercept)-Eira barbara                    -1.9932 0.7638 -3.5520 -1.9924
(Intercept)-Herpailurus yagouaroundi        -4.2024 1.4525 -7.0260 -4.2391
(Intercept)-Leopardus pardalis              -1.5218 0.3763 -2.3396 -1.4991
(Intercept)-Leopardus tigrinus              -6.9693 1.1498 -9.3426 -6.9447
(Intercept)-Leopardus wiedii                -3.2074 1.0814 -5.2825 -3.2580
(Intercept)-Mazama americana                -2.6124 0.5827 -3.9639 -2.5256
(Intercept)-Mazama nemorivaga               -5.6219 0.9905 -7.7977 -5.5421
(Intercept)-Myoprocta pratti                -4.7108 0.8809 -6.6393 -4.6596
(Intercept)-Myrmecophaga tridactyla         -1.7860 0.9374 -3.5221 -1.8715
(Intercept)-Nasua nasua                     -1.9922 0.9136 -3.8450 -1.9965
(Intercept)-Panthera onca                   -1.9128 0.6800 -3.3508 -1.8930
(Intercept)-Pecari tajacu                    0.4194 0.7086 -1.0669  0.4871
(Intercept)-Priodontes maximus              -1.6390 1.1842 -3.7822 -1.7003
(Intercept)-Procyon cancrivorous            -4.9635 1.4579 -7.8134 -4.9880
(Intercept)-Puma concolor                   -2.9256 0.6688 -4.4458 -2.8550
(Intercept)-Puma yagouaroundi               -5.0322 1.6389 -8.0014 -5.1341
(Intercept)-Tamandua tetradactyla           -3.3869 1.1391 -5.6753 -3.4186
(Intercept)-Tapirus terrestris              -0.3091 0.4660 -1.2972 -0.2661
(Intercept)-Tayassu pecari                  -1.5618 0.5460 -2.7912 -1.5187
scale(elev)-Atelocynus microtis             -0.6795 0.4607 -1.6793 -0.6568
scale(elev)-Coendou prehensilis              0.2016 0.6051 -0.9566  0.1761
scale(elev)-Cuniculus paca                  -0.1456 0.2690 -0.6826 -0.1429
scale(elev)-Dasyprocta fuliginosa            0.8599 0.4573 -0.1191  0.8688
scale(elev)-Dasypus sp.                     -0.2859 0.2911 -0.8687 -0.2831
scale(elev)-Didelphis marsupialis            0.1060 0.4523 -0.7272  0.0949
scale(elev)-Eira barbara                    -0.9852 0.4773 -2.0155 -0.9490
scale(elev)-Herpailurus yagouaroundi        -0.4436 0.6313 -1.7201 -0.4200
scale(elev)-Leopardus pardalis               0.3041 0.3235 -0.3392  0.3005
scale(elev)-Leopardus tigrinus              -0.2349 0.6160 -1.4665 -0.2163
scale(elev)-Leopardus wiedii                -0.1222 0.4773 -1.0382 -0.1275
scale(elev)-Mazama americana                 1.2261 0.3805  0.5285  1.2061
scale(elev)-Mazama nemorivaga                0.5441 0.3884 -0.1647  0.5344
scale(elev)-Myoprocta pratti                 0.0385 0.4762 -0.8527  0.0181
scale(elev)-Myrmecophaga tridactyla         -0.4303 0.4509 -1.3868 -0.4068
scale(elev)-Nasua nasua                     -0.2548 0.4632 -1.1907 -0.2454
scale(elev)-Panthera onca                    0.1494 0.4052 -0.6333  0.1401
scale(elev)-Pecari tajacu                    0.6037 0.3436 -0.0837  0.5997
scale(elev)-Priodontes maximus              -0.2231 0.5187 -1.2885 -0.2090
scale(elev)-Procyon cancrivorous            -0.1784 0.5066 -1.1401 -0.1878
scale(elev)-Puma concolor                   -0.3868 0.3827 -1.1346 -0.3870
scale(elev)-Puma yagouaroundi               -0.3696 0.5411 -1.4343 -0.3815
scale(elev)-Tamandua tetradactyla           -0.0505 0.5550 -1.2236 -0.0482
scale(elev)-Tapirus terrestris               0.2339 0.2977 -0.3225  0.2220
scale(elev)-Tayassu pecari                  -0.4401 0.3188 -1.0760 -0.4358
scale(border_dist)-Atelocynus microtis      -0.1040 0.4601 -0.9701 -0.1124
scale(border_dist)-Coendou prehensilis      -0.6611 0.9902 -2.7998 -0.5842
scale(border_dist)-Cuniculus paca           -1.2495 0.3798 -2.0256 -1.2391
scale(border_dist)-Dasyprocta fuliginosa    -0.5637 0.7374 -1.8886 -0.6360
scale(border_dist)-Dasypus sp.              -4.0608 0.8583 -5.8878 -4.0247
scale(border_dist)-Didelphis marsupialis     1.2427 0.5531  0.1607  1.2166
scale(border_dist)-Eira barbara             -1.6165 0.6216 -2.9293 -1.5937
scale(border_dist)-Herpailurus yagouaroundi  0.0974 0.9840 -1.8233  0.0702
scale(border_dist)-Leopardus pardalis        1.1100 0.3651  0.4592  1.0941
scale(border_dist)-Leopardus tigrinus       -0.0452 0.8988 -1.8211 -0.0628
scale(border_dist)-Leopardus wiedii          0.4129 0.5755 -0.6960  0.3801
scale(border_dist)-Mazama americana          0.7601 0.3915 -0.0159  0.7666
scale(border_dist)-Mazama nemorivaga        -3.1633 0.8865 -5.0055 -3.1105
scale(border_dist)-Myoprocta pratti         -0.2613 0.6597 -1.5601 -0.2578
scale(border_dist)-Myrmecophaga tridactyla  -0.3333 0.5332 -1.5022 -0.3107
scale(border_dist)-Nasua nasua              -0.7209 0.6225 -2.0352 -0.7064
scale(border_dist)-Panthera onca             0.5991 0.5032 -0.3344  0.5795
scale(border_dist)-Pecari tajacu            -0.4068 0.4605 -1.2674 -0.4185
scale(border_dist)-Priodontes maximus        0.2239 0.7064 -1.1391  0.1874
scale(border_dist)-Procyon cancrivorous     -1.7161 0.9004 -3.6384 -1.6456
scale(border_dist)-Puma concolor            -0.0506 0.4652 -0.9202 -0.0691
scale(border_dist)-Puma yagouaroundi        -1.3313 0.8941 -3.2487 -1.2827
scale(border_dist)-Tamandua tetradactyla    -0.1388 0.7359 -1.5393 -0.1511
scale(border_dist)-Tapirus terrestris       -0.2207 0.3598 -0.9291 -0.2216
scale(border_dist)-Tayassu pecari            0.2105 0.3731 -0.4775  0.1936
scale(FLII)-Atelocynus microtis              0.2162 0.4043 -0.5393  0.1876
scale(FLII)-Coendou prehensilis              0.5833 0.8316 -0.9115  0.5376
scale(FLII)-Cuniculus paca                  -0.3398 0.2082 -0.7608 -0.3320
scale(FLII)-Dasyprocta fuliginosa           -1.0294 0.3186 -1.6546 -1.0241
scale(FLII)-Dasypus sp.                      0.5703 0.6903 -0.6869  0.5438
scale(FLII)-Didelphis marsupialis            1.1550 0.5256  0.2919  1.0959
scale(FLII)-Eira barbara                    -0.4560 0.4952 -1.4884 -0.4208
scale(FLII)-Herpailurus yagouaroundi         0.4745 0.7454 -0.9289  0.4466
scale(FLII)-Leopardus pardalis               1.1397 0.3771  0.4929  1.1076
scale(FLII)-Leopardus tigrinus               0.7411 0.8301 -0.6506  0.6506
scale(FLII)-Leopardus wiedii                 0.9992 0.6638 -0.1296  0.9479
scale(FLII)-Mazama americana                 0.4242 0.2891 -0.0963  0.4028
scale(FLII)-Mazama nemorivaga                0.8028 0.7990 -0.5547  0.7441
scale(FLII)-Myoprocta pratti                 0.7582 0.5834 -0.2595  0.7023
scale(FLII)-Myrmecophaga tridactyla          0.3715 0.4431 -0.4923  0.3575
scale(FLII)-Nasua nasua                     -0.5512 0.5117 -1.7382 -0.4974
scale(FLII)-Panthera onca                    1.4655 0.6269  0.3808  1.4191
scale(FLII)-Pecari tajacu                    0.5874 0.2874  0.0654  0.5701
scale(FLII)-Priodontes maximus               0.2011 0.5526 -0.9146  0.2042
scale(FLII)-Procyon cancrivorous             0.8248 0.8854 -0.8008  0.7860
scale(FLII)-Puma concolor                    1.3410 0.6105  0.3110  1.2911
scale(FLII)-Puma yagouaroundi                0.6554 0.8666 -0.9249  0.5852
scale(FLII)-Tamandua tetradactyla           -0.6893 0.5788 -1.9848 -0.6390
scale(FLII)-Tapirus terrestris               1.0290 0.3384  0.4447  0.9974
scale(FLII)-Tayassu pecari                   1.7948 0.5805  0.8339  1.7362
                                              97.5%   Rhat  ESS
(Intercept)-Atelocynus microtis              0.1630 3.4943   77
(Intercept)-Coendou prehensilis             -1.6162 1.3372  115
(Intercept)-Cuniculus paca                   0.9270 4.4260    7
(Intercept)-Dasyprocta fuliginosa            0.6509 5.2340    7
(Intercept)-Dasypus sp.                     -3.6647 1.0256  366
(Intercept)-Didelphis marsupialis           -2.3853 1.9953   61
(Intercept)-Eira barbara                    -0.5427 1.3295  201
(Intercept)-Herpailurus yagouaroundi        -1.1533 1.2503  116
(Intercept)-Leopardus pardalis              -0.8511 1.0177  474
(Intercept)-Leopardus tigrinus              -4.8671 1.0263  314
(Intercept)-Leopardus wiedii                -0.9601 1.4396  152
(Intercept)-Mazama americana                -1.7187 1.4813  202
(Intercept)-Mazama nemorivaga               -3.9009 1.1896  396
(Intercept)-Myoprocta pratti                -3.1450 2.0065  118
(Intercept)-Myrmecophaga tridactyla          0.2959 1.0422  187
(Intercept)-Nasua nasua                     -0.1634 1.6135  130
(Intercept)-Panthera onca                   -0.6570 1.4455  197
(Intercept)-Pecari tajacu                    1.6999 4.2850    9
(Intercept)-Priodontes maximus               0.8206 1.2089  105
(Intercept)-Procyon cancrivorous            -2.0570 1.2222  101
(Intercept)-Puma concolor                   -1.7845 1.8071   57
(Intercept)-Puma yagouaroundi               -1.5904 1.0926  108
(Intercept)-Tamandua tetradactyla           -1.0210 1.1089  145
(Intercept)-Tapirus terrestris               0.5633 3.7843    8
(Intercept)-Tayassu pecari                  -0.6291 2.9047   15
scale(elev)-Atelocynus microtis              0.1768 1.0828 1135
scale(elev)-Coendou prehensilis              1.4312 1.0489 1254
scale(elev)-Cuniculus paca                   0.3539 1.0833  321
scale(elev)-Dasyprocta fuliginosa            1.7494 1.4850  102
scale(elev)-Dasypus sp.                      0.2698 1.0106 2104
scale(elev)-Didelphis marsupialis            1.0594 1.0311 1230
scale(elev)-Eira barbara                    -0.1166 1.0133 1030
scale(elev)-Herpailurus yagouaroundi         0.7641 1.0631 1340
scale(elev)-Leopardus pardalis               0.9502 1.0323 1124
scale(elev)-Leopardus tigrinus               0.9440 1.0123 1918
scale(elev)-Leopardus wiedii                 0.8115 1.0315 1815
scale(elev)-Mazama americana                 2.0137 1.1087  603
scale(elev)-Mazama nemorivaga                1.3553 1.0433 1111
scale(elev)-Myoprocta pratti                 1.0252 1.1106  462
scale(elev)-Myrmecophaga tridactyla          0.4233 1.0858 1038
scale(elev)-Nasua nasua                      0.6588 1.0470 1403
scale(elev)-Panthera onca                    0.9567 1.0482  880
scale(elev)-Pecari tajacu                    1.2818 1.2664  239
scale(elev)-Priodontes maximus               0.7682 1.0431 1348
scale(elev)-Procyon cancrivorous             0.8498 1.0006  825
scale(elev)-Puma concolor                    0.3668 1.1162  604
scale(elev)-Puma yagouaroundi                0.7427 1.0175 1006
scale(elev)-Tamandua tetradactyla            1.0254 1.0131 1412
scale(elev)-Tapirus terrestris               0.8563 1.2274  432
scale(elev)-Tayassu pecari                   0.1469 1.2317  257
scale(border_dist)-Atelocynus microtis       0.8099 1.1152  715
scale(border_dist)-Coendou prehensilis       1.1582 1.0117  318
scale(border_dist)-Cuniculus paca           -0.5146 1.5205   54
scale(border_dist)-Dasyprocta fuliginosa     1.0604 2.0477   33
scale(border_dist)-Dasypus sp.              -2.5126 1.0497  383
scale(border_dist)-Didelphis marsupialis     2.3785 1.2809  133
scale(border_dist)-Eira barbara             -0.4721 1.0355  665
scale(border_dist)-Herpailurus yagouaroundi  2.1216 1.0534  369
scale(border_dist)-Leopardus pardalis        1.8740 1.0002 1813
scale(border_dist)-Leopardus tigrinus        1.7372 1.0238  861
scale(border_dist)-Leopardus wiedii          1.5794 1.0611  694
scale(border_dist)-Mazama americana          1.5115 1.0763  325
scale(border_dist)-Mazama nemorivaga        -1.5573 1.1611  413
scale(border_dist)-Myoprocta pratti          1.0760 1.2949  137
scale(border_dist)-Myrmecophaga tridactyla   0.6404 1.0078  767
scale(border_dist)-Nasua nasua               0.4468 1.2030  335
scale(border_dist)-Panthera onca             1.6249 1.0939  848
scale(border_dist)-Pecari tajacu             0.5492 1.6689   39
scale(border_dist)-Priodontes maximus        1.6807 1.1737  384
scale(border_dist)-Procyon cancrivorous     -0.1425 1.0299  419
scale(border_dist)-Puma concolor             0.9134 1.2439  126
scale(border_dist)-Puma yagouaroundi         0.3016 1.0288  406
scale(border_dist)-Tamandua tetradactyla     1.3166 1.0413  741
scale(border_dist)-Tapirus terrestris        0.4983 1.2439  146
scale(border_dist)-Tayassu pecari            0.9672 1.3912  123
scale(FLII)-Atelocynus microtis              1.0503 1.1043  802
scale(FLII)-Coendou prehensilis              2.3472 1.0023  775
scale(FLII)-Cuniculus paca                   0.0495 1.0041 1565
scale(FLII)-Dasyprocta fuliginosa           -0.4221 1.0228  428
scale(FLII)-Dasypus sp.                      1.9990 1.0017 1259
scale(FLII)-Didelphis marsupialis            2.3503 1.0243  783
scale(FLII)-Eira barbara                     0.4375 1.0083  729
scale(FLII)-Herpailurus yagouaroundi         2.0757 1.0171  807
scale(FLII)-Leopardus pardalis               1.9778 1.0104 1202
scale(FLII)-Leopardus tigrinus               2.5449 1.0108  885
scale(FLII)-Leopardus wiedii                 2.3924 1.0072  685
scale(FLII)-Mazama americana                 1.0472 1.1076  672
scale(FLII)-Mazama nemorivaga                2.5582 1.0187  986
scale(FLII)-Myoprocta pratti                 2.0611 1.0038 1105
scale(FLII)-Myrmecophaga tridactyla          1.2919 1.0151 1022
scale(FLII)-Nasua nasua                      0.2977 1.1048  331
scale(FLII)-Panthera onca                    2.8124 1.0077  773
scale(FLII)-Pecari tajacu                    1.2021 1.0521  423
scale(FLII)-Priodontes maximus               1.3048 1.0220  703
scale(FLII)-Procyon cancrivorous             2.6493 1.0065  561
scale(FLII)-Puma concolor                    2.6847 1.0248  797
scale(FLII)-Puma yagouaroundi                2.5027 1.0045  679
scale(FLII)-Tamandua tetradactyla            0.3112 1.0159  562
scale(FLII)-Tapirus terrestris               1.7777 1.0181  896
scale(FLII)-Tayassu pecari                   3.1378 1.0019  964

Detection (logit scale): 
                                          Mean     SD    2.5%     50%   97.5%
(Intercept)-Atelocynus microtis        -3.0201 0.8089 -4.7187 -2.8665 -1.7653
(Intercept)-Coendou prehensilis        -3.7274 1.2270 -6.3725 -3.6024 -1.6674
(Intercept)-Cuniculus paca             -0.7034 0.0588 -0.8191 -0.7035 -0.5914
(Intercept)-Dasyprocta fuliginosa      -0.4486 0.0685 -0.5890 -0.4480 -0.3146
(Intercept)-Dasypus sp.                -1.5297 0.1277 -1.7825 -1.5291 -1.2816
(Intercept)-Didelphis marsupialis      -2.2225 0.3585 -2.9726 -2.2016 -1.5703
(Intercept)-Eira barbara               -3.0061 0.2288 -3.4562 -3.0043 -2.5640
(Intercept)-Herpailurus yagouaroundi   -3.6227 0.8534 -5.2392 -3.6211 -2.0292
(Intercept)-Leopardus pardalis         -1.7960 0.1609 -2.1121 -1.7923 -1.4875
(Intercept)-Leopardus tigrinus         -1.2974 0.9703 -3.3709 -1.2317  0.4554
(Intercept)-Leopardus wiedii           -3.2309 0.6724 -4.5121 -3.2208 -1.9812
(Intercept)-Mazama americana           -1.1923 0.1145 -1.4242 -1.1893 -0.9702
(Intercept)-Mazama nemorivaga          -2.0094 0.1825 -2.3611 -2.0102 -1.6529
(Intercept)-Myoprocta pratti           -1.3795 0.2433 -1.8825 -1.3702 -0.9100
(Intercept)-Myrmecophaga tridactyla    -3.0316 0.4590 -3.9425 -3.0033 -2.1732
(Intercept)-Nasua nasua                -3.0199 0.2886 -3.5948 -3.0179 -2.4751
(Intercept)-Panthera onca              -2.4552 0.2404 -2.9369 -2.4485 -1.9991
(Intercept)-Pecari tajacu              -1.2419 0.0737 -1.3870 -1.2414 -1.1001
(Intercept)-Priodontes maximus         -3.5773 0.4800 -4.4927 -3.5784 -2.6447
(Intercept)-Procyon cancrivorous       -3.6308 0.8324 -5.3690 -3.5858 -2.1633
(Intercept)-Puma concolor              -1.7334 0.2056 -2.1609 -1.7266 -1.3528
(Intercept)-Puma yagouaroundi          -3.7902 0.9880 -5.7765 -3.7237 -2.0481
(Intercept)-Tamandua tetradactyla      -3.4662 0.5853 -4.6464 -3.4358 -2.3856
(Intercept)-Tapirus terrestris         -1.1063 0.0773 -1.2595 -1.1051 -0.9584
(Intercept)-Tayassu pecari             -1.6780 0.1457 -1.9735 -1.6751 -1.3953
scale(effort)-Atelocynus microtis       0.3415 0.1828 -0.0082  0.3371  0.7149
scale(effort)-Coendou prehensilis       0.3957 0.2171 -0.0271  0.3906  0.8344
scale(effort)-Cuniculus paca            0.3395 0.0689  0.2064  0.3390  0.4786
scale(effort)-Dasyprocta fuliginosa     0.4121 0.0878  0.2441  0.4110  0.5868
scale(effort)-Dasypus sp.               0.4412 0.1350  0.1980  0.4351  0.7216
scale(effort)-Didelphis marsupialis     0.5480 0.1950  0.2086  0.5335  0.9748
scale(effort)-Eira barbara              0.5033 0.1823  0.1712  0.4921  0.8985
scale(effort)-Herpailurus yagouaroundi  0.2705 0.2097 -0.1570  0.2766  0.6637
scale(effort)-Leopardus pardalis        0.4150 0.1203  0.1924  0.4105  0.6712
scale(effort)-Leopardus tigrinus        0.4375 0.2185  0.0179  0.4318  0.9001
scale(effort)-Leopardus wiedii          0.4474 0.2043  0.0718  0.4374  0.8667
scale(effort)-Mazama americana          0.5246 0.1384  0.2784  0.5159  0.8180
scale(effort)-Mazama nemorivaga         0.2845 0.1489 -0.0016  0.2822  0.5880
scale(effort)-Myoprocta pratti          0.3487 0.1752  0.0044  0.3503  0.6903
scale(effort)-Myrmecophaga tridactyla   0.3854 0.1720  0.0649  0.3807  0.7334
scale(effort)-Nasua nasua               0.5000 0.1883  0.1735  0.4866  0.9150
scale(effort)-Panthera onca             0.3084 0.1399  0.0511  0.3042  0.5927
scale(effort)-Pecari tajacu             0.3535 0.0792  0.2010  0.3510  0.5135
scale(effort)-Priodontes maximus        0.2302 0.1712 -0.1158  0.2342  0.5585
scale(effort)-Procyon cancrivorous      0.4267 0.2114  0.0164  0.4232  0.8711
scale(effort)-Puma concolor             0.3910 0.1459  0.1150  0.3864  0.6919
scale(effort)-Puma yagouaroundi         0.3999 0.2152 -0.0226  0.3906  0.8372
scale(effort)-Tamandua tetradactyla     0.2815 0.1951 -0.1207  0.2870  0.6659
scale(effort)-Tapirus terrestris        0.3858 0.0818  0.2291  0.3842  0.5543
scale(effort)-Tayassu pecari            0.5449 0.1330  0.3087  0.5333  0.8264
                                         Rhat  ESS
(Intercept)-Atelocynus microtis        3.4286   34
(Intercept)-Coendou prehensilis        1.3341  106
(Intercept)-Cuniculus paca             1.0034 3000
(Intercept)-Dasyprocta fuliginosa      1.0095 3000
(Intercept)-Dasypus sp.                1.0040 2972
(Intercept)-Didelphis marsupialis      1.0229  830
(Intercept)-Eira barbara               1.0168  671
(Intercept)-Herpailurus yagouaroundi   1.2309  164
(Intercept)-Leopardus pardalis         1.0017 1934
(Intercept)-Leopardus tigrinus         1.0050 1211
(Intercept)-Leopardus wiedii           1.1798  196
(Intercept)-Mazama americana           1.0041 3000
(Intercept)-Mazama nemorivaga          1.0025 2389
(Intercept)-Myoprocta pratti           1.0024 2808
(Intercept)-Myrmecophaga tridactyla    1.0396  232
(Intercept)-Nasua nasua                1.1009  432
(Intercept)-Panthera onca              1.0169  490
(Intercept)-Pecari tajacu              1.0515 2207
(Intercept)-Priodontes maximus         1.1538  164
(Intercept)-Procyon cancrivorous       1.0350  101
(Intercept)-Puma concolor              1.0227 1394
(Intercept)-Puma yagouaroundi          1.0985  102
(Intercept)-Tamandua tetradactyla      1.0793  203
(Intercept)-Tapirus terrestris         1.0000 2853
(Intercept)-Tayassu pecari             1.0003 1810
scale(effort)-Atelocynus microtis      1.0074 2629
scale(effort)-Coendou prehensilis      1.0015 3000
scale(effort)-Cuniculus paca           1.0008 3000
scale(effort)-Dasyprocta fuliginosa    1.0001 2970
scale(effort)-Dasypus sp.              1.0035 3000
scale(effort)-Didelphis marsupialis    1.0007 2736
scale(effort)-Eira barbara             1.0032 2334
scale(effort)-Herpailurus yagouaroundi 1.0007 3000
scale(effort)-Leopardus pardalis       1.0005 2795
scale(effort)-Leopardus tigrinus       1.0028 2801
scale(effort)-Leopardus wiedii         1.0044 2529
scale(effort)-Mazama americana         1.0018 2760
scale(effort)-Mazama nemorivaga        1.0018 3000
scale(effort)-Myoprocta pratti         1.0069 3000
scale(effort)-Myrmecophaga tridactyla  0.9997 2372
scale(effort)-Nasua nasua              1.0027 2352
scale(effort)-Panthera onca            1.0071 2795
scale(effort)-Pecari tajacu            1.0007 3000
scale(effort)-Priodontes maximus       1.0016 2658
scale(effort)-Procyon cancrivorous     1.0025 3000
scale(effort)-Puma concolor            1.0014 3000
scale(effort)-Puma yagouaroundi        1.0032 3000
scale(effort)-Tamandua tetradactyla    1.0007 3000
scale(effort)-Tapirus terrestris       1.0012 2848
scale(effort)-Tayassu pecari           1.0016 3000

----------------------------------------
    Spatial Covariance
----------------------------------------
         Mean     SD   2.5%     50%   97.5%   Rhat  ESS
phi-1  8.9235 9.0901 0.0000  6.4319 26.4508 4.0224    8
phi-2 13.2789 7.9866 0.4498 13.0502 26.7775 1.0019 3000
phi-3  4.7494 8.0864 0.0000  0.0000 26.0644 9.1416    3
phi-4 13.5261 7.9445 0.6362 13.4657 26.8778 1.0001 2677
phi-5 13.6591 8.0601 0.6792 13.3846 26.8695 1.0007 3000
Code
summary(out.sp, level = 'both')

Call:
sfMsPGOcc(occ.formula = ~scale(elev) + scale(border_dist) + scale(FLII), 
    det.formula = ~scale(effort), data = data_list, cov.model = "exponential", 
    NNGP = TRUE, n.neighbors = 15, n.factors = 5, n.batch = 600, 
    batch.length = 25, n.omp.threads = 6, n.report = 1000, n.burn = 5000, 
    n.thin = 10, n.chains = 3)

Samples per Chain: 15000
Burn-in: 5000
Thinning Rate: 10
Number of Chains: 3
Total Posterior Samples: 3000
Run Time (min): 30.0153

----------------------------------------
    Community Level
----------------------------------------
Occurrence Means (logit scale): 
                      Mean     SD    2.5%     50%   97.5%   Rhat ESS
(Intercept)        -2.7611 0.4954 -3.7565 -2.7599 -1.7782 1.0807 277
scale(elev)        -0.0372 0.1944 -0.4196 -0.0382  0.3414 1.1540 850
scale(border_dist) -0.4589 0.3319 -1.1346 -0.4491  0.1698 1.0799 895
scale(FLII)         0.5170 0.2460  0.0567  0.5049  1.0430 1.0122 885

Occurrence Variances (logit scale): 
                     Mean     SD   2.5%    50%   97.5%   Rhat  ESS
(Intercept)        5.1766 2.1408 2.2988 4.7363 10.7314 1.1570  346
scale(elev)        0.4771 0.2324 0.1718 0.4252  1.0559 1.0385 1145
scale(border_dist) 2.0565 0.8866 0.8445 1.8913  4.2675 1.0374  438
scale(FLII)        0.9112 0.4685 0.3170 0.8046  2.0406 1.0008  496

Detection Means (logit scale): 
                 Mean     SD    2.5%     50%   97.5%   Rhat  ESS
(Intercept)   -2.2605 0.2700 -2.8176 -2.2461 -1.7495 1.1593  329
scale(effort)  0.3955 0.0653  0.2704  0.3935  0.5316 1.0033 2257

Detection Variances (logit scale): 
                Mean     SD   2.5%    50%  97.5%   Rhat  ESS
(Intercept)   1.5039 0.6865 0.6342 1.3541 3.2076 1.2673  106
scale(effort) 0.0454 0.0261 0.0160 0.0392 0.1110 1.0008 2566

----------------------------------------
    Species Level
----------------------------------------
Occurrence (logit scale): 
                                               Mean     SD    2.5%     50%
(Intercept)-Atelocynus microtis             -2.5173 1.0719 -3.9337 -2.7944
(Intercept)-Coendou prehensilis             -5.4512 1.7234 -8.7614 -5.5698
(Intercept)-Cuniculus paca                   0.1873 0.4071 -0.5579  0.1865
(Intercept)-Dasyprocta fuliginosa           -1.3041 1.1359 -3.5264 -1.2099
(Intercept)-Dasypus sp.                     -5.2100 0.9039 -7.1616 -5.1349
(Intercept)-Didelphis marsupialis           -3.9183 0.8388 -5.6388 -3.8827
(Intercept)-Eira barbara                    -1.9932 0.7638 -3.5520 -1.9924
(Intercept)-Herpailurus yagouaroundi        -4.2024 1.4525 -7.0260 -4.2391
(Intercept)-Leopardus pardalis              -1.5218 0.3763 -2.3396 -1.4991
(Intercept)-Leopardus tigrinus              -6.9693 1.1498 -9.3426 -6.9447
(Intercept)-Leopardus wiedii                -3.2074 1.0814 -5.2825 -3.2580
(Intercept)-Mazama americana                -2.6124 0.5827 -3.9639 -2.5256
(Intercept)-Mazama nemorivaga               -5.6219 0.9905 -7.7977 -5.5421
(Intercept)-Myoprocta pratti                -4.7108 0.8809 -6.6393 -4.6596
(Intercept)-Myrmecophaga tridactyla         -1.7860 0.9374 -3.5221 -1.8715
(Intercept)-Nasua nasua                     -1.9922 0.9136 -3.8450 -1.9965
(Intercept)-Panthera onca                   -1.9128 0.6800 -3.3508 -1.8930
(Intercept)-Pecari tajacu                    0.4194 0.7086 -1.0669  0.4871
(Intercept)-Priodontes maximus              -1.6390 1.1842 -3.7822 -1.7003
(Intercept)-Procyon cancrivorous            -4.9635 1.4579 -7.8134 -4.9880
(Intercept)-Puma concolor                   -2.9256 0.6688 -4.4458 -2.8550
(Intercept)-Puma yagouaroundi               -5.0322 1.6389 -8.0014 -5.1341
(Intercept)-Tamandua tetradactyla           -3.3869 1.1391 -5.6753 -3.4186
(Intercept)-Tapirus terrestris              -0.3091 0.4660 -1.2972 -0.2661
(Intercept)-Tayassu pecari                  -1.5618 0.5460 -2.7912 -1.5187
scale(elev)-Atelocynus microtis             -0.6795 0.4607 -1.6793 -0.6568
scale(elev)-Coendou prehensilis              0.2016 0.6051 -0.9566  0.1761
scale(elev)-Cuniculus paca                  -0.1456 0.2690 -0.6826 -0.1429
scale(elev)-Dasyprocta fuliginosa            0.8599 0.4573 -0.1191  0.8688
scale(elev)-Dasypus sp.                     -0.2859 0.2911 -0.8687 -0.2831
scale(elev)-Didelphis marsupialis            0.1060 0.4523 -0.7272  0.0949
scale(elev)-Eira barbara                    -0.9852 0.4773 -2.0155 -0.9490
scale(elev)-Herpailurus yagouaroundi        -0.4436 0.6313 -1.7201 -0.4200
scale(elev)-Leopardus pardalis               0.3041 0.3235 -0.3392  0.3005
scale(elev)-Leopardus tigrinus              -0.2349 0.6160 -1.4665 -0.2163
scale(elev)-Leopardus wiedii                -0.1222 0.4773 -1.0382 -0.1275
scale(elev)-Mazama americana                 1.2261 0.3805  0.5285  1.2061
scale(elev)-Mazama nemorivaga                0.5441 0.3884 -0.1647  0.5344
scale(elev)-Myoprocta pratti                 0.0385 0.4762 -0.8527  0.0181
scale(elev)-Myrmecophaga tridactyla         -0.4303 0.4509 -1.3868 -0.4068
scale(elev)-Nasua nasua                     -0.2548 0.4632 -1.1907 -0.2454
scale(elev)-Panthera onca                    0.1494 0.4052 -0.6333  0.1401
scale(elev)-Pecari tajacu                    0.6037 0.3436 -0.0837  0.5997
scale(elev)-Priodontes maximus              -0.2231 0.5187 -1.2885 -0.2090
scale(elev)-Procyon cancrivorous            -0.1784 0.5066 -1.1401 -0.1878
scale(elev)-Puma concolor                   -0.3868 0.3827 -1.1346 -0.3870
scale(elev)-Puma yagouaroundi               -0.3696 0.5411 -1.4343 -0.3815
scale(elev)-Tamandua tetradactyla           -0.0505 0.5550 -1.2236 -0.0482
scale(elev)-Tapirus terrestris               0.2339 0.2977 -0.3225  0.2220
scale(elev)-Tayassu pecari                  -0.4401 0.3188 -1.0760 -0.4358
scale(border_dist)-Atelocynus microtis      -0.1040 0.4601 -0.9701 -0.1124
scale(border_dist)-Coendou prehensilis      -0.6611 0.9902 -2.7998 -0.5842
scale(border_dist)-Cuniculus paca           -1.2495 0.3798 -2.0256 -1.2391
scale(border_dist)-Dasyprocta fuliginosa    -0.5637 0.7374 -1.8886 -0.6360
scale(border_dist)-Dasypus sp.              -4.0608 0.8583 -5.8878 -4.0247
scale(border_dist)-Didelphis marsupialis     1.2427 0.5531  0.1607  1.2166
scale(border_dist)-Eira barbara             -1.6165 0.6216 -2.9293 -1.5937
scale(border_dist)-Herpailurus yagouaroundi  0.0974 0.9840 -1.8233  0.0702
scale(border_dist)-Leopardus pardalis        1.1100 0.3651  0.4592  1.0941
scale(border_dist)-Leopardus tigrinus       -0.0452 0.8988 -1.8211 -0.0628
scale(border_dist)-Leopardus wiedii          0.4129 0.5755 -0.6960  0.3801
scale(border_dist)-Mazama americana          0.7601 0.3915 -0.0159  0.7666
scale(border_dist)-Mazama nemorivaga        -3.1633 0.8865 -5.0055 -3.1105
scale(border_dist)-Myoprocta pratti         -0.2613 0.6597 -1.5601 -0.2578
scale(border_dist)-Myrmecophaga tridactyla  -0.3333 0.5332 -1.5022 -0.3107
scale(border_dist)-Nasua nasua              -0.7209 0.6225 -2.0352 -0.7064
scale(border_dist)-Panthera onca             0.5991 0.5032 -0.3344  0.5795
scale(border_dist)-Pecari tajacu            -0.4068 0.4605 -1.2674 -0.4185
scale(border_dist)-Priodontes maximus        0.2239 0.7064 -1.1391  0.1874
scale(border_dist)-Procyon cancrivorous     -1.7161 0.9004 -3.6384 -1.6456
scale(border_dist)-Puma concolor            -0.0506 0.4652 -0.9202 -0.0691
scale(border_dist)-Puma yagouaroundi        -1.3313 0.8941 -3.2487 -1.2827
scale(border_dist)-Tamandua tetradactyla    -0.1388 0.7359 -1.5393 -0.1511
scale(border_dist)-Tapirus terrestris       -0.2207 0.3598 -0.9291 -0.2216
scale(border_dist)-Tayassu pecari            0.2105 0.3731 -0.4775  0.1936
scale(FLII)-Atelocynus microtis              0.2162 0.4043 -0.5393  0.1876
scale(FLII)-Coendou prehensilis              0.5833 0.8316 -0.9115  0.5376
scale(FLII)-Cuniculus paca                  -0.3398 0.2082 -0.7608 -0.3320
scale(FLII)-Dasyprocta fuliginosa           -1.0294 0.3186 -1.6546 -1.0241
scale(FLII)-Dasypus sp.                      0.5703 0.6903 -0.6869  0.5438
scale(FLII)-Didelphis marsupialis            1.1550 0.5256  0.2919  1.0959
scale(FLII)-Eira barbara                    -0.4560 0.4952 -1.4884 -0.4208
scale(FLII)-Herpailurus yagouaroundi         0.4745 0.7454 -0.9289  0.4466
scale(FLII)-Leopardus pardalis               1.1397 0.3771  0.4929  1.1076
scale(FLII)-Leopardus tigrinus               0.7411 0.8301 -0.6506  0.6506
scale(FLII)-Leopardus wiedii                 0.9992 0.6638 -0.1296  0.9479
scale(FLII)-Mazama americana                 0.4242 0.2891 -0.0963  0.4028
scale(FLII)-Mazama nemorivaga                0.8028 0.7990 -0.5547  0.7441
scale(FLII)-Myoprocta pratti                 0.7582 0.5834 -0.2595  0.7023
scale(FLII)-Myrmecophaga tridactyla          0.3715 0.4431 -0.4923  0.3575
scale(FLII)-Nasua nasua                     -0.5512 0.5117 -1.7382 -0.4974
scale(FLII)-Panthera onca                    1.4655 0.6269  0.3808  1.4191
scale(FLII)-Pecari tajacu                    0.5874 0.2874  0.0654  0.5701
scale(FLII)-Priodontes maximus               0.2011 0.5526 -0.9146  0.2042
scale(FLII)-Procyon cancrivorous             0.8248 0.8854 -0.8008  0.7860
scale(FLII)-Puma concolor                    1.3410 0.6105  0.3110  1.2911
scale(FLII)-Puma yagouaroundi                0.6554 0.8666 -0.9249  0.5852
scale(FLII)-Tamandua tetradactyla           -0.6893 0.5788 -1.9848 -0.6390
scale(FLII)-Tapirus terrestris               1.0290 0.3384  0.4447  0.9974
scale(FLII)-Tayassu pecari                   1.7948 0.5805  0.8339  1.7362
                                              97.5%   Rhat  ESS
(Intercept)-Atelocynus microtis              0.1630 3.4943   77
(Intercept)-Coendou prehensilis             -1.6162 1.3372  115
(Intercept)-Cuniculus paca                   0.9270 4.4260    7
(Intercept)-Dasyprocta fuliginosa            0.6509 5.2340    7
(Intercept)-Dasypus sp.                     -3.6647 1.0256  366
(Intercept)-Didelphis marsupialis           -2.3853 1.9953   61
(Intercept)-Eira barbara                    -0.5427 1.3295  201
(Intercept)-Herpailurus yagouaroundi        -1.1533 1.2503  116
(Intercept)-Leopardus pardalis              -0.8511 1.0177  474
(Intercept)-Leopardus tigrinus              -4.8671 1.0263  314
(Intercept)-Leopardus wiedii                -0.9601 1.4396  152
(Intercept)-Mazama americana                -1.7187 1.4813  202
(Intercept)-Mazama nemorivaga               -3.9009 1.1896  396
(Intercept)-Myoprocta pratti                -3.1450 2.0065  118
(Intercept)-Myrmecophaga tridactyla          0.2959 1.0422  187
(Intercept)-Nasua nasua                     -0.1634 1.6135  130
(Intercept)-Panthera onca                   -0.6570 1.4455  197
(Intercept)-Pecari tajacu                    1.6999 4.2850    9
(Intercept)-Priodontes maximus               0.8206 1.2089  105
(Intercept)-Procyon cancrivorous            -2.0570 1.2222  101
(Intercept)-Puma concolor                   -1.7845 1.8071   57
(Intercept)-Puma yagouaroundi               -1.5904 1.0926  108
(Intercept)-Tamandua tetradactyla           -1.0210 1.1089  145
(Intercept)-Tapirus terrestris               0.5633 3.7843    8
(Intercept)-Tayassu pecari                  -0.6291 2.9047   15
scale(elev)-Atelocynus microtis              0.1768 1.0828 1135
scale(elev)-Coendou prehensilis              1.4312 1.0489 1254
scale(elev)-Cuniculus paca                   0.3539 1.0833  321
scale(elev)-Dasyprocta fuliginosa            1.7494 1.4850  102
scale(elev)-Dasypus sp.                      0.2698 1.0106 2104
scale(elev)-Didelphis marsupialis            1.0594 1.0311 1230
scale(elev)-Eira barbara                    -0.1166 1.0133 1030
scale(elev)-Herpailurus yagouaroundi         0.7641 1.0631 1340
scale(elev)-Leopardus pardalis               0.9502 1.0323 1124
scale(elev)-Leopardus tigrinus               0.9440 1.0123 1918
scale(elev)-Leopardus wiedii                 0.8115 1.0315 1815
scale(elev)-Mazama americana                 2.0137 1.1087  603
scale(elev)-Mazama nemorivaga                1.3553 1.0433 1111
scale(elev)-Myoprocta pratti                 1.0252 1.1106  462
scale(elev)-Myrmecophaga tridactyla          0.4233 1.0858 1038
scale(elev)-Nasua nasua                      0.6588 1.0470 1403
scale(elev)-Panthera onca                    0.9567 1.0482  880
scale(elev)-Pecari tajacu                    1.2818 1.2664  239
scale(elev)-Priodontes maximus               0.7682 1.0431 1348
scale(elev)-Procyon cancrivorous             0.8498 1.0006  825
scale(elev)-Puma concolor                    0.3668 1.1162  604
scale(elev)-Puma yagouaroundi                0.7427 1.0175 1006
scale(elev)-Tamandua tetradactyla            1.0254 1.0131 1412
scale(elev)-Tapirus terrestris               0.8563 1.2274  432
scale(elev)-Tayassu pecari                   0.1469 1.2317  257
scale(border_dist)-Atelocynus microtis       0.8099 1.1152  715
scale(border_dist)-Coendou prehensilis       1.1582 1.0117  318
scale(border_dist)-Cuniculus paca           -0.5146 1.5205   54
scale(border_dist)-Dasyprocta fuliginosa     1.0604 2.0477   33
scale(border_dist)-Dasypus sp.              -2.5126 1.0497  383
scale(border_dist)-Didelphis marsupialis     2.3785 1.2809  133
scale(border_dist)-Eira barbara             -0.4721 1.0355  665
scale(border_dist)-Herpailurus yagouaroundi  2.1216 1.0534  369
scale(border_dist)-Leopardus pardalis        1.8740 1.0002 1813
scale(border_dist)-Leopardus tigrinus        1.7372 1.0238  861
scale(border_dist)-Leopardus wiedii          1.5794 1.0611  694
scale(border_dist)-Mazama americana          1.5115 1.0763  325
scale(border_dist)-Mazama nemorivaga        -1.5573 1.1611  413
scale(border_dist)-Myoprocta pratti          1.0760 1.2949  137
scale(border_dist)-Myrmecophaga tridactyla   0.6404 1.0078  767
scale(border_dist)-Nasua nasua               0.4468 1.2030  335
scale(border_dist)-Panthera onca             1.6249 1.0939  848
scale(border_dist)-Pecari tajacu             0.5492 1.6689   39
scale(border_dist)-Priodontes maximus        1.6807 1.1737  384
scale(border_dist)-Procyon cancrivorous     -0.1425 1.0299  419
scale(border_dist)-Puma concolor             0.9134 1.2439  126
scale(border_dist)-Puma yagouaroundi         0.3016 1.0288  406
scale(border_dist)-Tamandua tetradactyla     1.3166 1.0413  741
scale(border_dist)-Tapirus terrestris        0.4983 1.2439  146
scale(border_dist)-Tayassu pecari            0.9672 1.3912  123
scale(FLII)-Atelocynus microtis              1.0503 1.1043  802
scale(FLII)-Coendou prehensilis              2.3472 1.0023  775
scale(FLII)-Cuniculus paca                   0.0495 1.0041 1565
scale(FLII)-Dasyprocta fuliginosa           -0.4221 1.0228  428
scale(FLII)-Dasypus sp.                      1.9990 1.0017 1259
scale(FLII)-Didelphis marsupialis            2.3503 1.0243  783
scale(FLII)-Eira barbara                     0.4375 1.0083  729
scale(FLII)-Herpailurus yagouaroundi         2.0757 1.0171  807
scale(FLII)-Leopardus pardalis               1.9778 1.0104 1202
scale(FLII)-Leopardus tigrinus               2.5449 1.0108  885
scale(FLII)-Leopardus wiedii                 2.3924 1.0072  685
scale(FLII)-Mazama americana                 1.0472 1.1076  672
scale(FLII)-Mazama nemorivaga                2.5582 1.0187  986
scale(FLII)-Myoprocta pratti                 2.0611 1.0038 1105
scale(FLII)-Myrmecophaga tridactyla          1.2919 1.0151 1022
scale(FLII)-Nasua nasua                      0.2977 1.1048  331
scale(FLII)-Panthera onca                    2.8124 1.0077  773
scale(FLII)-Pecari tajacu                    1.2021 1.0521  423
scale(FLII)-Priodontes maximus               1.3048 1.0220  703
scale(FLII)-Procyon cancrivorous             2.6493 1.0065  561
scale(FLII)-Puma concolor                    2.6847 1.0248  797
scale(FLII)-Puma yagouaroundi                2.5027 1.0045  679
scale(FLII)-Tamandua tetradactyla            0.3112 1.0159  562
scale(FLII)-Tapirus terrestris               1.7777 1.0181  896
scale(FLII)-Tayassu pecari                   3.1378 1.0019  964

Detection (logit scale): 
                                          Mean     SD    2.5%     50%   97.5%
(Intercept)-Atelocynus microtis        -3.0201 0.8089 -4.7187 -2.8665 -1.7653
(Intercept)-Coendou prehensilis        -3.7274 1.2270 -6.3725 -3.6024 -1.6674
(Intercept)-Cuniculus paca             -0.7034 0.0588 -0.8191 -0.7035 -0.5914
(Intercept)-Dasyprocta fuliginosa      -0.4486 0.0685 -0.5890 -0.4480 -0.3146
(Intercept)-Dasypus sp.                -1.5297 0.1277 -1.7825 -1.5291 -1.2816
(Intercept)-Didelphis marsupialis      -2.2225 0.3585 -2.9726 -2.2016 -1.5703
(Intercept)-Eira barbara               -3.0061 0.2288 -3.4562 -3.0043 -2.5640
(Intercept)-Herpailurus yagouaroundi   -3.6227 0.8534 -5.2392 -3.6211 -2.0292
(Intercept)-Leopardus pardalis         -1.7960 0.1609 -2.1121 -1.7923 -1.4875
(Intercept)-Leopardus tigrinus         -1.2974 0.9703 -3.3709 -1.2317  0.4554
(Intercept)-Leopardus wiedii           -3.2309 0.6724 -4.5121 -3.2208 -1.9812
(Intercept)-Mazama americana           -1.1923 0.1145 -1.4242 -1.1893 -0.9702
(Intercept)-Mazama nemorivaga          -2.0094 0.1825 -2.3611 -2.0102 -1.6529
(Intercept)-Myoprocta pratti           -1.3795 0.2433 -1.8825 -1.3702 -0.9100
(Intercept)-Myrmecophaga tridactyla    -3.0316 0.4590 -3.9425 -3.0033 -2.1732
(Intercept)-Nasua nasua                -3.0199 0.2886 -3.5948 -3.0179 -2.4751
(Intercept)-Panthera onca              -2.4552 0.2404 -2.9369 -2.4485 -1.9991
(Intercept)-Pecari tajacu              -1.2419 0.0737 -1.3870 -1.2414 -1.1001
(Intercept)-Priodontes maximus         -3.5773 0.4800 -4.4927 -3.5784 -2.6447
(Intercept)-Procyon cancrivorous       -3.6308 0.8324 -5.3690 -3.5858 -2.1633
(Intercept)-Puma concolor              -1.7334 0.2056 -2.1609 -1.7266 -1.3528
(Intercept)-Puma yagouaroundi          -3.7902 0.9880 -5.7765 -3.7237 -2.0481
(Intercept)-Tamandua tetradactyla      -3.4662 0.5853 -4.6464 -3.4358 -2.3856
(Intercept)-Tapirus terrestris         -1.1063 0.0773 -1.2595 -1.1051 -0.9584
(Intercept)-Tayassu pecari             -1.6780 0.1457 -1.9735 -1.6751 -1.3953
scale(effort)-Atelocynus microtis       0.3415 0.1828 -0.0082  0.3371  0.7149
scale(effort)-Coendou prehensilis       0.3957 0.2171 -0.0271  0.3906  0.8344
scale(effort)-Cuniculus paca            0.3395 0.0689  0.2064  0.3390  0.4786
scale(effort)-Dasyprocta fuliginosa     0.4121 0.0878  0.2441  0.4110  0.5868
scale(effort)-Dasypus sp.               0.4412 0.1350  0.1980  0.4351  0.7216
scale(effort)-Didelphis marsupialis     0.5480 0.1950  0.2086  0.5335  0.9748
scale(effort)-Eira barbara              0.5033 0.1823  0.1712  0.4921  0.8985
scale(effort)-Herpailurus yagouaroundi  0.2705 0.2097 -0.1570  0.2766  0.6637
scale(effort)-Leopardus pardalis        0.4150 0.1203  0.1924  0.4105  0.6712
scale(effort)-Leopardus tigrinus        0.4375 0.2185  0.0179  0.4318  0.9001
scale(effort)-Leopardus wiedii          0.4474 0.2043  0.0718  0.4374  0.8667
scale(effort)-Mazama americana          0.5246 0.1384  0.2784  0.5159  0.8180
scale(effort)-Mazama nemorivaga         0.2845 0.1489 -0.0016  0.2822  0.5880
scale(effort)-Myoprocta pratti          0.3487 0.1752  0.0044  0.3503  0.6903
scale(effort)-Myrmecophaga tridactyla   0.3854 0.1720  0.0649  0.3807  0.7334
scale(effort)-Nasua nasua               0.5000 0.1883  0.1735  0.4866  0.9150
scale(effort)-Panthera onca             0.3084 0.1399  0.0511  0.3042  0.5927
scale(effort)-Pecari tajacu             0.3535 0.0792  0.2010  0.3510  0.5135
scale(effort)-Priodontes maximus        0.2302 0.1712 -0.1158  0.2342  0.5585
scale(effort)-Procyon cancrivorous      0.4267 0.2114  0.0164  0.4232  0.8711
scale(effort)-Puma concolor             0.3910 0.1459  0.1150  0.3864  0.6919
scale(effort)-Puma yagouaroundi         0.3999 0.2152 -0.0226  0.3906  0.8372
scale(effort)-Tamandua tetradactyla     0.2815 0.1951 -0.1207  0.2870  0.6659
scale(effort)-Tapirus terrestris        0.3858 0.0818  0.2291  0.3842  0.5543
scale(effort)-Tayassu pecari            0.5449 0.1330  0.3087  0.5333  0.8264
                                         Rhat  ESS
(Intercept)-Atelocynus microtis        3.4286   34
(Intercept)-Coendou prehensilis        1.3341  106
(Intercept)-Cuniculus paca             1.0034 3000
(Intercept)-Dasyprocta fuliginosa      1.0095 3000
(Intercept)-Dasypus sp.                1.0040 2972
(Intercept)-Didelphis marsupialis      1.0229  830
(Intercept)-Eira barbara               1.0168  671
(Intercept)-Herpailurus yagouaroundi   1.2309  164
(Intercept)-Leopardus pardalis         1.0017 1934
(Intercept)-Leopardus tigrinus         1.0050 1211
(Intercept)-Leopardus wiedii           1.1798  196
(Intercept)-Mazama americana           1.0041 3000
(Intercept)-Mazama nemorivaga          1.0025 2389
(Intercept)-Myoprocta pratti           1.0024 2808
(Intercept)-Myrmecophaga tridactyla    1.0396  232
(Intercept)-Nasua nasua                1.1009  432
(Intercept)-Panthera onca              1.0169  490
(Intercept)-Pecari tajacu              1.0515 2207
(Intercept)-Priodontes maximus         1.1538  164
(Intercept)-Procyon cancrivorous       1.0350  101
(Intercept)-Puma concolor              1.0227 1394
(Intercept)-Puma yagouaroundi          1.0985  102
(Intercept)-Tamandua tetradactyla      1.0793  203
(Intercept)-Tapirus terrestris         1.0000 2853
(Intercept)-Tayassu pecari             1.0003 1810
scale(effort)-Atelocynus microtis      1.0074 2629
scale(effort)-Coendou prehensilis      1.0015 3000
scale(effort)-Cuniculus paca           1.0008 3000
scale(effort)-Dasyprocta fuliginosa    1.0001 2970
scale(effort)-Dasypus sp.              1.0035 3000
scale(effort)-Didelphis marsupialis    1.0007 2736
scale(effort)-Eira barbara             1.0032 2334
scale(effort)-Herpailurus yagouaroundi 1.0007 3000
scale(effort)-Leopardus pardalis       1.0005 2795
scale(effort)-Leopardus tigrinus       1.0028 2801
scale(effort)-Leopardus wiedii         1.0044 2529
scale(effort)-Mazama americana         1.0018 2760
scale(effort)-Mazama nemorivaga        1.0018 3000
scale(effort)-Myoprocta pratti         1.0069 3000
scale(effort)-Myrmecophaga tridactyla  0.9997 2372
scale(effort)-Nasua nasua              1.0027 2352
scale(effort)-Panthera onca            1.0071 2795
scale(effort)-Pecari tajacu            1.0007 3000
scale(effort)-Priodontes maximus       1.0016 2658
scale(effort)-Procyon cancrivorous     1.0025 3000
scale(effort)-Puma concolor            1.0014 3000
scale(effort)-Puma yagouaroundi        1.0032 3000
scale(effort)-Tamandua tetradactyla    1.0007 3000
scale(effort)-Tapirus terrestris       1.0012 2848
scale(effort)-Tayassu pecari           1.0016 3000

----------------------------------------
    Spatial Covariance
----------------------------------------
         Mean     SD   2.5%     50%   97.5%   Rhat  ESS
phi-1  8.9235 9.0901 0.0000  6.4319 26.4508 4.0224    8
phi-2 13.2789 7.9866 0.4498 13.0502 26.7775 1.0019 3000
phi-3  4.7494 8.0864 0.0000  0.0000 26.0644 9.1416    3
phi-4 13.5261 7.9445 0.6362 13.4657 26.8778 1.0001 2677
phi-5 13.6591 8.0601 0.6792 13.3846 26.8695 1.0007 3000

Bayesian p-values can be inspected to check for lack of fit (overall or by species). Lack of fit at significance level = 0.05 is indicated by Bayesian p-values below 0.025 or greater than 0.975. The overall Bayesian p-value (Bpvalue) indicates no problems with lack of fit. Likewise, species-level Bayesian p-values (Bpvalue_species) indicate no lack of fit for any species.

Gelman and Rubin’s convergence diagnostic: Approximate convergence is diagnosed when the upper limit is close to 1.

Convergence is diagnosed when the chains have ‘forgotten’ their initial values, and the output from all chains is indistinguishable. The gelman.diag diagnostic is applied to a single variable from the chain. It is based a comparison of within-chain and between-chain variances, and is similar to a classical analysis of variance.

Values substantially above 1 indicate lack of convergence.

Model Diagnostics

Code
#| eval: true
#| echo: true
#| code-fold: true
#| warning: false
#| message: false
# Extract posterior draws for later use
posterior1 <- as.array(out.sp)

# Trace plots to check chain mixing. Extract posterior samples and bind in a single matrix.
POSTERIOR.MATRIX <- cbind(out.sp$alpha.comm.samples, 
                          out.sp$beta.comm.samples,  
                          out.sp$alpha.samples, 
                          out.sp$beta.samples)

# Matrix output is all chains combined, split into 3 chains.
CHAIN.1 <- as.mcmc(POSTERIOR.MATRIX[1:1000,])
CHAIN.2 <- as.mcmc(POSTERIOR.MATRIX[1001:2000,])
CHAIN.3 <- as.mcmc(POSTERIOR.MATRIX[2001:3000,])
# CHAIN.4 <- as.mcmc(POSTERIOR.MATRIX[8001:10000,])

# Bind four chains as coda mcmc.list object.
POSTERIOR.CHAINS <- mcmc.list(CHAIN.1, CHAIN.2, CHAIN.3)#, CHAIN.4)

# Create an empty folder.
# dir.create ("Beetle_plots")

# Plot chain mixing of each parameter to a multi-panel plot and save to the new folder. ART 5 mins

######################################
#### SAVE Diagnostics at PDF
# MCMCtrace(POSTERIOR.CHAINS, params = "all", Rhat = TRUE, n.eff = TRUE)#, pdf = TRUE, filename = "Beetle_240909_traceplots.pdf", wd = "Beetle_plots")



#mcmc_trace(fit.commu, parms = c("beta.ranef.cont.border_dist.mean"))

#posterior2 <- extract(fit.commu, inc_warmup = TRUE, permuted = FALSE)

#color_scheme_set("mix-blue-pink")
#p <- mcmc_trace(posterior1,  pars = c("mu", "tau"), n_warmup = 300,
#                facet_args = list(nrow = 2, labeller = label_parsed))
#p + facet_text(size = 15)



#outMCMC <- fit.commu #Convert output to MCMC object
#diagnostics chains 

# all as pdf
# MCMCtrace(outMCMC)

# MCMCtrace(outMCMC, params = c("alpha0"), type = 'trace', Rhat = TRUE, n.eff = TRUE)

# MCMCtrace(outMCMC, params = c("beta0"), type = 'trace', Rhat = TRUE, n.eff = TRUE)

# MCMCtrace(outMCMC, params = c("beta.ranef.cont.border_dist"), type = 'trace', Rhat = TRUE, n.eff = TRUE)

# MCMCtrace(out.sp, params = c("beta.ranef.cont.border_dist.mean"), type = 'trace', pdf = F, Rhat = TRUE, n.eff = TRUE)

# MCMCtrace(outMCMC, params = c("beta.ranef.cont.elev"), type = 'trace', Rhat = TRUE, n.eff = TRUE)

# MCMCtrace(outMCMC, params = c("beta.ranef.cont.elev.mean"), type = 'trace', pdf = F, Rhat = TRUE, n.eff = TRUE)


### density
# MCMCtrace(outMCMC, params = c("Nspecies"), ISB = FALSE, pdf = F, exact = TRUE, post_zm = TRUE, type = 'density', Rhat = TRUE, n.eff = TRUE, ind = TRUE)

### density
# MCMCtrace(outMCMC, params = c("beta.ranef.cont.elev.mean"), ISB = FALSE, pdf = F, exact = TRUE, post_zm = TRUE, type = 'density', Rhat = TRUE, n.eff = TRUE, ind = TRUE)

### density
#MCMCtrace(outMCMC, params = c("beta.ranef.cont.border_dist.mean"), ISB = FALSE, pdf = F, exact = TRUE, post_zm = TRUE, type = 'density', Rhat = TRUE, n.eff = TRUE, ind = TRUE)

#coda::gelman.diag(outMCMC,  multivariate = FALSE, transform=FALSE)
                    
# coda::gelman.plot(outMCMC,  multivariate = FALSE)



# 
# mcmc_intervals(outMCMC, pars = c("Nspecies_in_AP[1]",
#                                  "Nspecies_in_AP[2]"),
#                point_est = "mean",
#                prob = 0.75, prob_outer = 0.95) +
#   ggtitle("Number of species") + 
#   scale_y_discrete(labels = c("Nspecies_in_AP[1]"=levels(sitecovs$in_AP)[1],
#              "Nspecies_in_AP[2]"=levels(sitecovs$in_AP)[2]))
# 
# #Continuous
# p <- mcmc_intervals(outMCMC, 
#                pars = c("beta.ranef.cont.border_dist.mean",
#                          #"beta.ranef.cont.elev.mean",
#                         "beta.ranef.categ.in_AP.mean[2]"))
# 
# # relabel parameters
# p + scale_y_discrete(
#   labels = c("beta.ranef.cont.border_dist.mean"="Dist_border",
#                          #"beta.ranef.cont.elev.mean"="Elevation",
#                         "beta.ranef.categ.in_AP.mean[2]"="in_AP")
# ) +
#   ggtitle("Treatment effect on all species")
# 

Posterior Summary (effect)

Community effects

Code
# Create simple plot summaries using MCMCvis package.
# Detection covariate effects --------- 
MCMCplot(out.sp$alpha.comm.samples, ref_ovl = TRUE, ci = c(50, 95))
# Occupancy community-level effects 
MCMCplot(out.sp$beta.comm.samples, ref_ovl = TRUE, ci = c(50, 95))
(a) Detection
(b) Occupancy
Figure 1: Community effects

Species effects

Code
# Occupancy species-level effects 
MCMCplot(out.sp$beta.samples[,26:50], ref_ovl = TRUE, ci = c(50, 95))

# Occupancy species-level effects 
MCMCplot(out.sp$beta.samples[,51:75], ref_ovl = TRUE, ci = c(50, 95))
(a) elevation
(b) distance border
Figure 2: Species effects

Species effects using bayesplot

Code
library(bayesplot)
#mcmc_areas(outMCMC, regex_pars = "Nspecies_in_AP")
# mcmc_areas(outMCMC, regex_pars = "Nspecies_in_AP")

mcmc_intervals(out.sp$beta.samples[,26:50] , point_est = "mean",
               prob = 0.75, prob_outer = 0.95) + 
  geom_vline(xintercept = 0, color = "red", linetype = "dashed", size = 0.5)
Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.

Code
mcmc_intervals(out.sp$beta.samples[,51:75] , point_est = "mean",
               prob = 0.75, prob_outer = 0.95) + 
  geom_vline(xintercept = 0, color = "red", linetype = "dashed", size = 0.5)

Code
mcmc_intervals(out.sp$beta.samples[,76:100] , point_est = "mean",
               prob = 0.75, prob_outer = 0.95) + 
  geom_vline(xintercept = 0, color = "red", linetype = "dashed", size = 0.5)

Prediction as graph

This prediction uses the non spatial model

Code
# 6. Prediction -----------------------------------------------------------
# Predict occupancy along a gradient of elev   
# Create a set of values across the range of observed elev values
elev.pred.vals <- seq(min(data_list$occ.covs$elev), 
                              max(data_list$occ.covs$elev), 
                              length.out = 100)
# Scale predicted values by mean and standard deviation used to fit the model
elev.pred.vals.scale <- (elev.pred.vals - mean(data_list$occ.covs$elev)) / 
                             sd(data_list$occ.covs$elev)
# Create a set of values across the range of observed elev values
border_dist.pred.vals <- seq(min(data_list$occ.covs$border_dist), 
                              max(data_list$occ.covs$border_dist), 
                              length.out = 100)
# Scale predicted values by mean and standard deviation used to fit the model
border_dist.pred.vals.scale <- (border_dist.pred.vals -
                                  mean(data_list$occ.covs$border_dist)) /
                                  sd(data_list$occ.covs$border_dist)

# Create a set of values across the range of observed FLII values
FLII.pred.vals <- seq(min(data_list$occ.covs$FLII), 
                              max(data_list$occ.covs$FLII), 
                              length.out = 100)
# Scale predicted values by mean and standard deviation used to fit the model
FLII.pred.vals.scale <- (FLII.pred.vals -
                                  mean(data_list$occ.covs$FLII)) /
                                  sd(data_list$occ.covs$FLII)

# Predict occupancy across elev  values at mean values of all other variables
pred.df1 <- as.matrix(data.frame(intercept = 1, 
                                 elev = elev.pred.vals.scale, 
                                 border_dist = 0,
                                 FLII=0))#, catchment = 0, density = 0, 
                         # slope = 0))
# Predict occupancy across elev  values at mean values of all other variables
pred.df2 <- as.matrix(data.frame(intercept = 1, elev = 0, 
                         border_dist = border_dist.pred.vals.scale,
                         FLII= 0 ))#, catchment = 0, density = 0, 
                         # slope = 0))
# Predict occupancy across elev  values at mean values of all other variables
pred.df3 <- as.matrix(data.frame(intercept = 1, elev = 0, 
                         border_dist = 0,
                         FLII= FLII.pred.vals.scale ))#, catchment = 0, density = 0, 
                         # slope = 0))

out.pred1 <- predict(out, pred.df1) # using non spatial
str(out.pred1)
List of 3
 $ psi.0.samples: num [1:1500, 1:25, 1:100] 0.0664 0.1062 0.1333 0.1399 0.1003 ...
 $ z.0.samples  : int [1:1500, 1:25, 1:100] 0 0 0 0 1 1 0 0 1 0 ...
 $ call         : language predict.msPGOcc(object = out, X.0 = pred.df1)
 - attr(*, "class")= chr "predict.msPGOcc"
Code
psi.0.quants <- apply(out.pred1$psi.0.samples, c(2, 3), quantile, 
                          prob = c(0.025, 0.5, 0.975))
sp.codes <- attr(data_list$y, "dimnames")[[1]]
psi.plot.dat <- data.frame(psi.med = c(t(psi.0.quants[2, , ])), 
                                 psi.low = c(t(psi.0.quants[1, , ])), 
                                 psi.high = c(t(psi.0.quants[3, , ])), 
                           elev = elev.pred.vals, 
                                 sp.codes = rep(sp.codes, 
                                                each = length(elev.pred.vals)))

ggplot(psi.plot.dat, aes(x = elev, y = psi.med)) + 
  geom_ribbon(aes(ymin = psi.low, ymax = psi.high), fill = 'grey70') +
  geom_line() + 
  facet_wrap(vars(sp.codes)) + 
  theme_bw() + 
  labs(x = 'elevation (m)', y = 'Occupancy Probability') 

Code
out.pred2 <- predict(out, pred.df2) # using non spatial
str(out.pred2)
List of 3
 $ psi.0.samples: num [1:1500, 1:25, 1:100] 0.0726 0.1975 0.15 0.1059 0.0332 ...
 $ z.0.samples  : int [1:1500, 1:25, 1:100] 0 1 0 0 0 0 0 0 0 0 ...
 $ call         : language predict.msPGOcc(object = out, X.0 = pred.df2)
 - attr(*, "class")= chr "predict.msPGOcc"
Code
psi.0.quants <- apply(out.pred2$psi.0.samples, c(2, 3), quantile, 
                          prob = c(0.025, 0.5, 0.975))
sp.codes <- attr(data_list$y, "dimnames")[[1]]
psi.plot.dat <- data.frame(psi.med = c(t(psi.0.quants[2, , ])), 
                                 psi.low = c(t(psi.0.quants[1, , ])), 
                                 psi.high = c(t(psi.0.quants[3, , ])), 
                           border_dist = border_dist.pred.vals, 
                                 sp.codes = rep(sp.codes, 
                                                each = length(border_dist.pred.vals)))

ggplot(psi.plot.dat, aes(x = border_dist/1000, y = psi.med)) + 
  geom_ribbon(aes(ymin = psi.low, ymax = psi.high), fill = 'grey70') +
  geom_line() + 
  facet_wrap(vars(sp.codes)) + 
  theme_bw() + 
  labs(x = 'border_dist (Km)', y = 'Occupancy Probability') 

Code
out.pred3 <- predict(out, pred.df3) # using non spatial
str(out.pred3)
List of 3
 $ psi.0.samples: num [1:1500, 1:25, 1:100] 0.0348 0.14116 0.09034 0.01307 0.00479 ...
 $ z.0.samples  : int [1:1500, 1:25, 1:100] 0 0 0 0 0 0 0 0 0 0 ...
 $ call         : language predict.msPGOcc(object = out, X.0 = pred.df3)
 - attr(*, "class")= chr "predict.msPGOcc"
Code
psi.0.quants <- apply(out.pred3$psi.0.samples, c(2, 3), quantile, 
                          prob = c(0.025, 0.5, 0.975))
sp.codes <- attr(data_list$y, "dimnames")[[1]]
psi.plot.dat <- data.frame(psi.med = c(t(psi.0.quants[2, , ])), 
                                 psi.low = c(t(psi.0.quants[1, , ])), 
                                 psi.high = c(t(psi.0.quants[3, , ])), 
                           FLII = FLII.pred.vals, 
                                 sp.codes = rep(sp.codes, 
                                                each = length(FLII.pred.vals)))

ggplot(psi.plot.dat, aes(x = FLII/1000, y = psi.med)) + 
  geom_ribbon(aes(ymin = psi.low, ymax = psi.high), fill = 'grey70') +
  geom_line() + 
  facet_wrap(vars(sp.codes)) + 
  theme_bw() + 
  labs(x = 'FLII', y = 'Occupancy Probability') 

Spatial Prediction in elevation_EC

Code
#aggregate  resolution to (factor = 3)
#transform coord data to UTM
elevation_UTM <- project(elevation_EC, "EPSG:10603")
elevation_17.aggregate <- aggregate(elevation_UTM, fact=10)
res(elevation_17.aggregate)
[1] 5984.313 5984.313
Code
# Convert the SpatRaster to a data frame with coordinates
df_coords <- as.data.frame(elevation_17.aggregate, xy = TRUE)
names(df_coords) <-c("X","Y","elev")

elev.pred <- (df_coords$elev - mean(data_list$occ.covs$elev)) / sd(data_list$occ.covs$elev)


############### Predict new data ############### 
# we predict at one covariate varing and others at mean value = 0
# intercept = 1, elev = 0, border_dist = 0, FLII= 0
# in that order
################ ################ ################ 
predict_data <- cbind(1, elev.pred, 0, 0)
colnames(predict_data) <- c("intercept",
                            "elev",
                            "border_dist",
                            "FLII" )

# X.0 <- cbind(1, 1, elev.pred)#, elev.pred^2)

# coords.0 <- as.matrix(hbefElev[, c('Easting', 'Northing')])
out.sp.ms.pred <- predict(out.sp, 
                          X.0=predict_data, 
                          df_coords[,1:2],
                          threads=4) #Warning :'threads' is not an argument
Warning in predict.sfMsPGOcc(out.sp, X.0 = predict_data, df_coords[, 1:2], :
'threads' is not an argument
----------------------------------------
    Prediction description
----------------------------------------
Spatial Factor NNGP Multispecies Occupancy model with Polya-Gamma latent
variable fit with 428 observations.

Number of covariates 4 (including intercept if specified).

Using the exponential spatial correlation model.

Using 15 nearest neighbors.
Using 5 latent spatial factors.

Number of MCMC samples 3000.

Predicting at 2618 non-sampled locations.


Source compiled with OpenMP support and model fit using 1 threads.
-------------------------------------------------
        Predicting
-------------------------------------------------
Location: 100 of 2618, 3.82%
Location: 200 of 2618, 7.64%
Location: 300 of 2618, 11.46%
Location: 400 of 2618, 15.28%
Location: 500 of 2618, 19.10%
Location: 600 of 2618, 22.92%
Location: 700 of 2618, 26.74%
Location: 800 of 2618, 30.56%
Location: 900 of 2618, 34.38%
Location: 1000 of 2618, 38.20%
Location: 1100 of 2618, 42.02%
Location: 1200 of 2618, 45.84%
Location: 1300 of 2618, 49.66%
Location: 1400 of 2618, 53.48%
Location: 1500 of 2618, 57.30%
Location: 1600 of 2618, 61.12%
Location: 1700 of 2618, 64.94%
Location: 1800 of 2618, 68.75%
Location: 1900 of 2618, 72.57%
Location: 2000 of 2618, 76.39%
Location: 2100 of 2618, 80.21%
Location: 2200 of 2618, 84.03%
Location: 2300 of 2618, 87.85%
Location: 2400 of 2618, 91.67%
Location: 2500 of 2618, 95.49%
Location: 2600 of 2618, 99.31%
Location: 2618 of 2618, 100.00%
Generating latent occupancy state
Code
# extract the array of interest= occupancy
predicted_array <- out.sp.ms.pred$psi.0.samples



dim(predicted_array)
[1] 3000   25 2618
Code
# df_plot <- NA # empty column
predicted_raster_list <- list() # rast(nrows = 52, ncols = 103,
                         #ext(elevation_17.aggregate), 
                         # crs = "EPSG:32717")

###################################################
# simpler way to make mean in the array by species
# array order: itera=3000, sp=22, pixels=5202
# we wan to keep index 2 and 3
library(plyr)
------------------------------------------------------------------------------
You have loaded plyr after dplyr - this is likely to cause problems.
If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
library(plyr); library(dplyr)
------------------------------------------------------------------------------

Adjuntando el paquete: 'plyr'
The following objects are masked from 'package:dplyr':

    arrange, count, desc, failwith, id, mutate, rename, summarise,
    summarize
The following object is masked from 'package:purrr':

    compact
The following object is masked from 'package:maps':

    ozone
Code
aresult = (plyr::aaply(predicted_array, # mean fo all iterations
                       c(2,3), 
                       mean))

##############################################
#### Loop to make rasters and put in a list
# array order: itera=3000, sp=22, pixels=5202
for (i in 1:dim(predicted_array)[2]){
 # Convert array slice to data frame
 # df_plot[i] <- as.data.frame(predicted_array[1, i, ]) # sp_i
 # df_plot[i] <- as.vector(aresult[i,]) # Extract sp
 # Producing an SDM for all (posterior mean)
 plot.dat <- data.frame(x = df_coords$X, 
                        y = df_coords$Y, 
                        psi.sp = as.vector(aresult[i,]))
 pred_rast <- rast(plot.dat, 
                   type = "xyz", 
                   crs = "EPSG:10603") # Replace EPSG:4326 with your CRS
 predicted_raster_list[[i]] <- pred_rast
}

# Convert the list to a SpatRaster stack
predictad_raster_stack <- rast(predicted_raster_list)
# put species names
names(predictad_raster_stack) <- selected.sp

plot(predictad_raster_stack)

Code
# get the mean
predicted_mean <- mean(predictad_raster_stack)

plot(predicted_mean, main="mean occupancy")

Code
mapview(predicted_mean) + mapview(AP_Yasuni_UTM_line)
Code
# 

Package Citation

Code
pkgs <- cite_packages(output = "paragraph", pkgs="Session", out.dir = ".")
# knitr::kable(pkgs)
pkgs

We used R v. 4.4.2 (R Core Team 2024) and the following R packages: abind v. 1.4.8 (Plate and Heiberger 2024), bayesplot v. 1.14.0 (Gabry et al. 2019; Gabry and Mahr 2025), beepr v. 2.0 (Bååth 2024), camtrapR v. 3.0.0 (Niedballa et al. 2016), coda v. 0.19.4.1 (Plummer et al. 2006), DT v. 0.34.0 (Xie et al. 2025), elevatr v. 0.99.0 (Hollister et al. 2023), maps v. 3.4.3 (Becker et al. 2025), mapview v. 2.11.4 (Appelhans et al. 2025), MCMCvis v. 0.16.3 (Youngflesh 2018), plyr v. 1.8.9 (Wickham 2011), sf v. 1.0.21 (Pebesma 2018; Pebesma and Bivand 2023), snow v. 0.4.4 (Tierney et al. 2021), snowfall v. 1.84.6.3 (Knaus 2023), spOccupancy v. 0.8.0 (Doser et al. 2022, 2024; Doser, Finley, and Banerjee 2023), terra v. 1.8.70 (Hijmans 2025), tictoc v. 1.2.1 (Izrailev 2024), tidyverse v. 2.0.0 (Wickham et al. 2019), tmap v. 4.2 (Tennekes 2018).

Sesion info

Code
print(sessionInfo(), locale = FALSE)
R version 4.4.2 (2024-10-31 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 10 x64 (build 19045)

Matrix products: internal


attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] plyr_1.8.9        abind_1.4-8       lubridate_1.9.4   forcats_1.0.0    
 [5] stringr_1.5.2     dplyr_1.1.4       purrr_1.1.0       readr_2.1.5      
 [9] tidyr_1.3.1       tibble_3.2.1      ggplot2_4.0.0     tidyverse_2.0.0  
[13] spOccupancy_0.8.0 camtrapR_3.0.0    snowfall_1.84-6.3 snow_0.4-4       
[17] beepr_2.0         coda_0.19-4.1     MCMCvis_0.16.3    tictoc_1.2.1     
[21] bayesplot_1.14.0  elevatr_0.99.0    terra_1.8-70      tmap_4.2         
[25] maps_3.4.3        mapview_2.11.4    sf_1.0-21         DT_0.34.0        
[29] readxl_1.4.3      grateful_0.3.0   

loaded via a namespace (and not attached):
  [1] RColorBrewer_1.1-3      tensorA_0.36.2.1        rstudioapi_0.17.1      
  [4] audio_0.1-11            jsonlite_2.0.0          wk_0.9.4               
  [7] magrittr_2.0.3          farver_2.1.2            nloptr_2.1.1           
 [10] rmarkdown_2.30          vctrs_0.6.5             minqa_1.2.8            
 [13] base64enc_0.1-3         RcppNumerical_0.6-0     progress_1.2.3         
 [16] htmltools_0.5.8.1       leafsync_0.1.0          distributional_0.5.0   
 [19] curl_7.0.0              raster_3.6-32           cellranger_1.1.0       
 [22] s2_1.1.9                sass_0.4.10             bslib_0.9.0            
 [25] slippymath_0.3.1        KernSmooth_2.23-24      htmlwidgets_1.6.4      
 [28] cachem_1.1.0            stars_0.6-8             uuid_1.2-1             
 [31] mime_0.13               lifecycle_1.0.4         iterators_1.0.14       
 [34] pkgconfig_2.0.3         cols4all_0.8-1          Matrix_1.7-1           
 [37] R6_2.6.1                fastmap_1.2.0           shiny_1.9.1            
 [40] digest_0.6.37           colorspace_2.1-1        leafem_0.2.4           
 [43] crosstalk_1.2.1         labeling_0.4.3          lwgeom_0.2-14          
 [46] progressr_0.15.0        spacesXYZ_1.6-0         timechange_0.3.0       
 [49] httr_1.4.7              mgcv_1.9-1              compiler_4.4.2         
 [52] microbenchmark_1.5.0    proxy_0.4-27            withr_3.0.2            
 [55] doParallel_1.0.17       backports_1.5.0         brew_1.0-10            
 [58] S7_0.2.0                DBI_1.2.3               logger_0.4.0           
 [61] MASS_7.3-61             maptiles_0.10.0         tmaptools_3.3          
 [64] leaflet_2.2.3           classInt_0.4-11         tools_4.4.2            
 [67] units_0.8-7             leaflegend_1.2.1        httpuv_1.6.16          
 [70] glue_1.8.0              satellite_1.0.5         nlme_3.1-166           
 [73] promises_1.3.3          grid_4.4.2              checkmate_2.3.2        
 [76] reshape2_1.4.4          generics_0.1.3          leaflet.providers_2.0.0
 [79] gtable_0.3.6            tzdb_0.4.0              shinyBS_0.61.1         
 [82] class_7.3-22            data.table_1.17.8       hms_1.1.3              
 [85] sp_2.2-0                RANN_2.6.2              foreach_1.5.2          
 [88] pillar_1.11.1           posterior_1.6.1         later_1.4.2            
 [91] splines_4.4.2           lattice_0.22-6          tidyselect_1.2.1       
 [94] knitr_1.50              svglite_2.1.3           stats4_4.4.2           
 [97] xfun_0.52               shinydashboard_0.7.3    leafpop_0.1.0          
[100] stringi_1.8.4           rematch_2.0.0           yaml_2.3.10            
[103] boot_1.3-31             evaluate_1.0.4          codetools_0.2-20       
[106] cli_3.6.5               RcppParallel_5.1.9      systemfonts_1.1.0      
[109] xtable_1.8-4            jquerylib_0.1.4         secr_5.1.0             
[112] dichromat_2.0-0.1       Rcpp_1.1.0              spAbundance_0.2.1      
[115] png_0.1-8               XML_3.99-0.18           parallel_4.4.2         
[118] prettyunits_1.2.0       lme4_1.1-35.5           mvtnorm_1.3-2          
[121] scales_1.4.0            e1071_1.7-16            crayon_1.5.3           
[124] rlang_1.1.6            

References

Appelhans, Tim, Florian Detsch, Christoph Reudenbach, and Stefan Woellauer. 2025. mapview: Interactive Viewing of Spatial Data in r. https://CRAN.R-project.org/package=mapview.
Bååth, Rasmus. 2024. beepr: Easily Play Notification Sounds on Any Platform. https://CRAN.R-project.org/package=beepr.
Becker, Richard A., Allan R. Wilks, Ray Brownrigg, Thomas P. Minka, and Alex Deckmyn. 2025. maps: Draw Geographical Maps. https://CRAN.R-project.org/package=maps.
Doser, Jeffrey W., Andrew O. Finley, and Sudipto Banerjee. 2023. “Joint Species Distribution Models with Imperfect Detection for High-Dimensional Spatial Data.” Ecology, e4137. https://doi.org/10.1002/ecy.4137.
Doser, Jeffrey W., Andrew O. Finley, Marc Kéry, and Elise F. Zipkin. 2022. spOccupancy: An r Package for Single-Species, Multi-Species, and Integrated Spatial Occupancy Models.” Methods in Ecology and Evolution 13: 1670–78. https://doi.org/10.1111/2041-210X.13897.
Doser, Jeffrey W., Andrew O. Finley, Sarah P. Saunders, Marc Kéry, Aaron S. Weed, and Elise F. Zipkin. 2024. “Modeling Complex Species-Environment Relationships Through Spatially-Varying Coefficient Occupancy Models.” Journal of Agricultural, Biological, and Environmental Statistics. https://doi.org/10.1007/s13253-023-00595-6.
Gabry, Jonah, and Tristan Mahr. 2025. bayesplot: Plotting for Bayesian Models.” https://mc-stan.org/bayesplot/.
Gabry, Jonah, Daniel Simpson, Aki Vehtari, Michael Betancourt, and Andrew Gelman. 2019. “Visualization in Bayesian Workflow.” J. R. Stat. Soc. A 182: 389–402. https://doi.org/10.1111/rssa.12378.
Hijmans, Robert J. 2025. terra: Spatial Data Analysis. https://CRAN.R-project.org/package=terra.
Hollister, Jeffrey, Tarak Shah, Jakub Nowosad, Alec L. Robitaille, Marcus W. Beck, and Mike Johnson. 2023. elevatr: Access Elevation Data from Various APIs. https://doi.org/10.5281/zenodo.8335450.
Izrailev, Sergei. 2024. tictoc: Functions for Timing r Scripts, as Well as Implementations of Stack and StackList Structures. https://CRAN.R-project.org/package=tictoc.
Knaus, Jochen. 2023. snowfall: Easier Cluster Computing (Based on snow). https://CRAN.R-project.org/package=snowfall.
Niedballa, Jürgen, Rahel Sollmann, Alexandre Courtiol, and Andreas Wilting. 2016. camtrapR: An r Package for Efficient Camera Trap Data Management.” Methods in Ecology and Evolution 7 (12): 1457–62. https://doi.org/10.1111/2041-210X.12600.
Pebesma, Edzer. 2018. Simple Features for R: Standardized Support for Spatial Vector Data.” The R Journal 10 (1): 439–46. https://doi.org/10.32614/RJ-2018-009.
Pebesma, Edzer, and Roger Bivand. 2023. Spatial Data Science: With applications in R. Chapman and Hall/CRC. https://doi.org/10.1201/9780429459016.
Plate, Tony, and Richard Heiberger. 2024. abind: Combine Multidimensional Arrays. https://CRAN.R-project.org/package=abind.
Plummer, Martyn, Nicky Best, Kate Cowles, and Karen Vines. 2006. CODA: Convergence Diagnosis and Output Analysis for MCMC.” R News 6 (1): 7–11. https://journal.r-project.org/archive/.
R Core Team. 2024. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.
Tennekes, Martijn. 2018. tmap: Thematic Maps in R.” Journal of Statistical Software 84 (6): 1–39. https://doi.org/10.18637/jss.v084.i06.
Tierney, Luke, A. J. Rossini, Na Li, and H. Sevcikova. 2021. snow: Simple Network of Workstations. https://CRAN.R-project.org/package=snow.
Wickham, Hadley. 2011. “The Split-Apply-Combine Strategy for Data Analysis.” Journal of Statistical Software 40 (1): 1–29. https://www.jstatsoft.org/v40/i01/.
Wickham, Hadley, Mara Averick, Jennifer Bryan, Winston Chang, Lucy D’Agostino McGowan, Romain François, Garrett Grolemund, et al. 2019. “Welcome to the tidyverse.” Journal of Open Source Software 4 (43): 1686. https://doi.org/10.21105/joss.01686.
Xie, Yihui, Joe Cheng, Xianying Tan, and Garrick Aden-Buie. 2025. DT: A Wrapper of the JavaScript Library DataTables. https://CRAN.R-project.org/package=DT.
Youngflesh, Casey. 2018. MCMCvis: Tools to Visualize, Manipulate, and Summarize MCMC Output.” Journal of Open Source Software 3 (24): 640. https://doi.org/10.21105/joss.00640.

Reuse

Citation

BibTeX citation:
@online{forero2025,
  author = {Forero, German and Wallace, Robert and Zapara-Rios, Galo and
    Isasi-Catalá, Emiliana and J. Lizcano, Diego},
  title = {Fitting a {Spatial} {Factor} {Multi-Species} {Occupancy}
    {Model}},
  date = {2025-08-16},
  url = {https://dlizcano.github.io/Occu_APs_all/blog/2025-10-15-analysis/},
  langid = {en}
}
For attribution, please cite this work as:
Forero, German, Robert Wallace, Galo Zapara-Rios, Emiliana Isasi-Catalá, and Diego J. Lizcano. 2025. “Fitting a Spatial Factor Multi-Species Occupancy Model.” August 16, 2025. https://dlizcano.github.io/Occu_APs_all/blog/2025-10-15-analysis/.